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ABSTRACT 

Problems  in  accurate  determination  of  seismic 
threshold  magnitude  are  analytically  tractable  for 
only  some  simple,  but  interesting,  cases.  This  report 
gives  these  analytic  results  and  proceeds  to  a simula- 
tion experiment  to  provide  insights  where  analytic 
predictions  are  not  possible.  Differences  in  direct, 
incremental,  and  cumulative  threshold  determinations; 
distinctions  between  true  and  observed  thresholds, 
effects  of  correlated  signals  and  noise;  feedback 
between  LR  and  P thresholds;  and  estimation  of  Ms_mb 
relationships  are  treated  in  this  simulation  experi- 
ment. A comparison  of  simulated  detection  results  is 
made  with  real  data  from  LASA  on  the  Kamchatka-Kurile 
region  and  good  agreement  is  obtained  between  predic- 
tions and  observations. 
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INTRODUCTION 


Determination  of  seismic  threshold  for  a station 
or  twork  of  stations  serves  to  evaluate  that  station 
or  letwork  in  the  context  of  seismic  detection.  Threshold 
iu /el  reflects  several  different  aspects  of  the 
station  or  network:  seismic  background  noise,  relative 

signal  strength,  seismometry,  recording  equipment, 
processing  techniques,  and  analysis  procedures.  If 
one  aspires  to  quantify  these  effects  or  to  evaluate 
station  and  network  detection  capability,  it  is  then 
necessary  to  make  an  accurate  determination  of  threshold. 
It  is  also  desirable  to  have  a technique  for 
predicting  thresholds  and  for  comparing  predictions  to 
observations . 

Determining  seismic  thresholds  has  several  facets 
in  practice.  An  empirical  threshold  can  be  determined 
directly  by  plotting  the  percentage  of  events  detected 
in  each  equal  magnitude  increment  against  the  known 
magnitudes.  This  is  the  direct  method;  it  gives  an 
unbiased  threshold  only  if  one  has  an  independent 
source  of  magnitude  information  whose  threshold  is  much 
lower  than  the  station  or  network  under  consideration 
and  whose  reported  magnitudes  have  a small  variance  and 
are  insignificantly  biased.  Since  an  independent  source 
of  such  nature  is  often  unavailable,  one  must  use  a 
method  dependent  upon  magnitude  information  from  the 
network  or  station  under  consideration;  this  involves 
plotting  the  lumber  of  events  having  a certain  phase 
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(e „ g . , LR)  detected  against  the  magnitudes  calculated 
from  that  phase  (e.g.,  Mc).  This  is  the  incremental 
method;  thresholds  are  found  by  estimating  the  true 
seismicity-magnitude  relation  from  the  asymptotic  part 
of  the  data  at  high  magnitudes  or  by  a maximum-likeli- 
hood method  (Kelley  and  Lacoss,  1969).  We  will  show 
that  the  use  of  observed,  rather  than  true  or  indepen- 
dent, magnitudes  distorts  the  simple  picture  presented 
by  the  direct  method  and  effects  the  threshold  deter- 
mination. A variation  of  the  incremental  method  is 
to  plot  the  cumulative  number  of  detections  (with 
magnitude  greater  than  or  equal  to  a given  observed 
magnitude)  against  the  observed  magnitudes;  this  is 
the  cumulative  method.  We  will  show  that  this  approach, 
taken  to  smooth  the  data  points,  heightens  the  threshold 
distortion  even  more. 

Another  somewhat  different  approach  is  used  in 
seismology  and  is  a type  of  direct  threshold  determina- 
tion. Here,  in  the  manner  of  the  direct  method,  the 
percentage  of  detections  of  surface  waves  in  a given 
body-wave  magnitude  increment  is  plotted  against  the 
body-wave  magnitude,  supplied  by  a source  which  pre- 
sumably detects  body  waves  much  more  frequently  than 
the  station  or  network  under  consideration  detects 
surface  vaves.  Thresholds  are  then  converted  from  the 

m,  to  M scale.  We  will  show  how  this  method  of  M 
b s s 

threshold  determination  is  so  fraught  with  pitfalls 
that  only  the  most  careful  attention  to  detail  can 
produce  reliable  thresholds. 
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The  matter  of  calculated  thresholds,  by  use  of 
program  NbTWORTII  (Wirth,  1970),  also  deserves  attention 
since  the  authors  have  noted  cases  of  unacceptable 
discrepancy  between  calculated  and  observed  thresholds. 
All  the  discrepancy  is  not  necessarily  due  to  faulty 
empirical  trreshold  determination,  and  we  must  recon- 
sider the  assumptions  in  the  common  scheme  for  calcu- 
lating station  or  network  thresholds. 

This  report  discusses  several  of  the  important 
aspects  of  threshold  determination,  calculated  and 
observed,  which  we  feel  have  not  been  given  exposure 
or  properly  interrelated;  we  feel  that  attention  to 
these  aspects,  along  with  those  presented  by  previous 
investigators,  will  result  in  unbiased  threshold 
determinations  which  are  suitable  for  use  in  comparisons 
between  stations  and  between  networks. 

We  must  make  some  simplifying  assumptions  about 
certain  phenomena  in  order  to  make  the  mathematics 
tractable;  they  are,  however,  realistic  and  commonly 
used.  Therefore,  we  allow  the  signal  amplitudes  from 
an  event  to  be  lognormally  distributed  (Freedman,  1967; 
von  Seggern,  1972).  The  same  applies  to  the  noise 
amplitudes  at  a station  (Geotech,  1966).  We  allow  the 
logarithm  of  the  number  of  events  which  occur  in  a 
specified  time  interval  to  be  proportional  to  magnitude 
(Gutenberg  and  Richter,  1954;  Scholz,  1968);  that  is, 
we  use  a linear  seismicity-amplitude  relation  through- 
out. We  also  assume  a linear  relationship  between  M 

, s 

and  m^  everywhere. 
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The  first  parts  of  this  report  establish  the 
concepts,  aiul  theoretical  derivations  are  given  when 
necessary.  A later  part  of  this  report  describes  the 
method  and  results  of  a simulation  experiment  which 
serves  to  verify  the  predictions  of  the  first  parts 
and  to  give  insights  where  mathematical  predictions 
are  not  attainable.  The  techniques  used  in  this 
simulation  are  essential  tools  lor  satisfactory 
evaluations  of  practical  networks.  The  last  part  of 
the  report  deals  with  refinements  in  calculating 
seismic  thresholds. 

To  aid  the  reader,  we  provide  next  for  rererence 
a list  of  notations  appearing  in  equations  throughout 
the  report  and  also  a glossary  of  important  recurrent 
terms  used  in  the  report.. 
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NOTATION 


M 


0 


II 


A 


m, 


M 


m. 


M 

M 


m 


mb 


a 

ms 


a 

s 


a 

ss 


mean  common  logarithm  of  the  noise  amplitude 
at  a station 

standard  deviation  of  the  common  logarithm 
of  noise  amplitude  at  a station 
calculated  amplitude  of  the  signal  at  a 
given  station 

observed  amplitude  of  the  signal  at  a given 
s tation 

true  body-wave  magnitude 
true  surface-wave  magnitude 
operational  body-wave  magnitude 
operational  surface-wave  magnitude 
observed  body-wave  magnitude 
observed  surface-wave  magnitude 
standard  deviation  of  operational  magnitudes 
about  the  true  magnitude 

standard  deviation  of  operational  body-wave 
magnitudes  about  true  magnitudes 
standard  deviation  of  operational  surface- 
wave  magnitudes  about  true  magnitudes 
standard  deviation  of  observed  magnitudes 
about  operational  magnitude 
standard  deviation  of  observed  body-wave 
magnitudes  about  operational  magnitudes 
standard  deviation  of  observed  surface-wave 
magnitudes  about  operational  magnitudes 


e 
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3 


:T 


a 

b 

7 


erf 

I’jjCm) 

I*!  Cm) 
Cm) 

P(x) 

1’Uly) 

“a 


f(x) 

f(x|y) 


number  or  events  occurring  with  magnitude  m 
number  of  events  detected  with  magnitude  m 
slope  of  the  base  1U  seismicity-magnitude 
re  la t ion 

slope  of  the  base  e seismicity-magnitude 
relation 

intercept  of  the  seismicity-magnitude 
re lat ion 
••  io6 

slope  of  the  M -m^  relation 
intercept  of  the  M - ni|  relation 
normal  probability  density  function 
unit  normal  probability  distribution 
funct  ion 
error  function 

direct  probability  of  detection  curve 
incremental  probability  of  detection  curve 
cumulative  probability  of  detection  curve 
probability  distribution  function 
conditional  probability  distribution  function 
percentage  point  for  the  unit  normal  dis- 
tribution 

binomial  probability  for  exactly  k of  n 
occur ranees 

probability  density  function 
conditional  probability  density  function 
correlation  coefficient 
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GLOSSARY 


Bias 

source  bias  - a difference  between  the 

operational  magnitude  and 
the  true  magnitude  due  to 
differences  in  fault  dimen- 
sions and  shape,  existing 
stress  levels,  close-in 
clastic  parameters,  etc., 
which  are  not  explicitly 
allowed  for  in  the  network's 
routine  determ ' nations  of 
magnitude . 

path-station  bias  - a difference  between  the 

observed  magnitude  and  the 
operational  magnitude  due  to 
differences  in  depth,  radia- 
tion pattern,  path  effects, 
station  characteristics,  etc., 
whose  effects  may  be  removed 
from  the  magnitude  estimates 
by  the  network's  routine 
analysis . 

threshold  bias  - a difference  between  observed 

and  true  threshold  magnitudes 
resulting  from  incremental 
and  cumulative  threshold 
determination  methods. 
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Magnitude 

t rue  magnitude 


operational  magnitude 


the  scale  which  most  nearly 
represents  the  energy  release 
of  a seismic  event;  there  is 
a one-to-one  mapping  between 
this  scale  and  e.g.,  "yield" 
for  explosions  or  "seismic 
moment"  for  earthquakes. 


tied  to  a particular  phase 
such  as  Ffm^)  or  LR(Ms) , 
this  is  the  magnitude  which 
an  ideal  network  (in  which 
a large  number  of  stations 
with  extremely  low  threshold 
densely  cover  the  globe), 
using  its  routine  highly- 
refined  magnitude  estimation 
techniques,  would  estimate 
for  a seismic  event. 


observed  magnitude  - the  magnitude  of  a seismic 

event  as  routinely  deter- 
mined by  a station  or  real 
network  (in  which  a limited 
number  of  stations  with 
varying  thresholds  sparsely 
cover  the  globe)  . 
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Correlation 


correlated  noise 


correlated  signals 


Relation 

seismicity- magnitude 
relation 


Ms-m^  relation 


Threshold 

threshold  magnitude 


noise  with  non-zero  correla- 
tion of  changes  in  long-term 
1 hour)  levels  among  sta- 
tions of  a network. 

correlated  differences  among 
stations  between  the  observed 
station  magnitudes  and  true 
magnitude.  These  differences 
are  related  to  source  bias. 

the  equation  expressing  the 
dependence  of  the  number  of 
events  occurring  within  a 
specified  time  on  magnitude, 
usually  of  the  form  log(num- 
ber)  = a(magnitude)  + B. 

the  equation  expressing  the 

dependence  of  surface-wave 

magnitude  on  body  wave 

magnitude  or  vice-versa, 

usually  of  the  foim 

M = am.  + b . 
s b 

a suitably  chosen  magnitude 
at  which  a certain  propor- 
tion, usually  50%  or  90%, 
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r 


direct  threshold 
dete  rmi nation 


r 

incremental  threshold 
determination 


of  the  evtents  of  interest 
are  detected  via  a phase 
such  a-'  I*  or  LR. 

a method  (?f  threshold  magni- 
tude determination  whereby 
the  percentage  of  events  for 
which  a phase  of  interest 
(e.g.,  I’  or  LR)  was  detected 
is  plotted  against  some 
magnitude  which  is  not 
calculated  from  those  phases 
themselves . 

a method  of  threshold  magni- 
tude determination  whereby 
the  number  of  event  detections 
via  a certain  phase  (e.g., 

P or  LR)  in  each  magnitude 
increment  is  plotted  against 
magnitudes  determined  usually, 
but  not  always,  from  these 
phases  themselves.  Thresh- 
olds are  determined  where 
plotted  points  begin  to 
deviate  from  the  straight- 
line  approximation  to  the 
seismicity-magnitude  rela- 
tion, which  is  determined 
by  higher-magnitude  events 
where  1001  detection  is 
achieved. 
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cumulative  threshold 
determination 


same  as  incremental  thresh 
old  determination,  c/cept 
that  cumulative  numbers  of 
detections  with  magnitudes 
greater  than  or  equal  to 
the  current  magnitude  are 
plotted. 
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CONCEPT  OF  SEISMIC  MAGNITUDE 


In  discussing  seismic  thresholds,  it  is  necessary 
to  view  seismic  magnitude  quite  precisely.  We  have  a 
need  for  three  conceptions  of  magnitud::  true,  opera- 

tional and  observed,  as  defined  below.  Such  a tricho- 
tomy lias  been  implied  for  explosion  magnitudes  by  von 
Seggern  (1973).  We  must  then  discuss  our  theoretical 
predictions  and  practical  applications  in  this  frame- 
work and  clearly  distinguish  at  all  times  which  magni- 
tude concept  we  are  using. 

The  question  of  the  ’'true"  magnitude  of  a seismic 
event  is  a difficult  one  for  which  no  completely 
satisfying  answer  is  available.  One  approach  to  this 
problem  is  to  regard  the  magnitude  which  maps  one-to- 
one  into  the  radiometric  or  design  yield  of  an  explo- 
sion as  the  fundamental,  or  as  we  shall  say  the  "true", 
magnitude;  for  earthquakes,  "true"  magnitude  could  be 
linked  to  seismic  moment  or  the  total  radiated  energy. 

Even  for  a perfect  unbiased  network,  the  "opera- 
tional" magnitudes,  such  as  m^ , M , or  M,  , as  computed 
by  standard  se ismological  procedures,  would  vary  for 
events  of  the  same  true  magnitude.  This  variation  we 
call  source  bias.  For  example,  difference  in  depths 
of  burial  or  in  properties  of  the  surrounding  media 
will  result  in  different  magnitudes  for  the  same  yield. 
Similarly,  different  fault  mechanisms,  depths  of  focus, 
and  states  of  stress  can  result  in  varying  operational 
magnitudes  for  otherwise  identical  earthquakes.  Finally, 
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we  define  our  real  estimates  of  magnitude  as  the 
"observed"  magnitudes;  these  will  vary  about  the 
operational  magnitude  due  to  radiation  patterns,  sta- 
tion noise  levels,  path  and  station  effects,  etc. 

The  true  magnitudes  are  seismo logically  unmeasurable 
and  therefore  unavailable  to  use.  Likewise  opera- 
tional magnitudes,  although  measurable  in  principle, 
arc  unavailable  to  use  because  of  the  limitations  of 
real  seismic  networks.  In  general,  the  larger  the 
c^ent,  the  more  closely  careful  analysis  will  cause 
the  observed  magnitude  to  approach  the  operational 
magnitude . 

I'or  presenting  results  it  will  be  desirable  to 
express  thresholds  in  terms  of  one  or  more  of  the  above 
three  magnitudes,  depending  on  the  frame  of  reference 
of  the  discussions.  for  example,  if  one  asks  for  the 
threshold  of  detection  of  a 2 kt  shot  in  granite,  the 
answer  is  different,  as  we  shall  show,  from  that  of 
an  event  with  operational  m^  equal  to  that  of  the 
expected  value  of  the  operational  m^  for  2 kt  in 
granite.  Of  course  fer  most  cases  of  interest  these 
differences  arc  not  great  (on  the  order  of  0.1  magni- 
tude unit)  . 
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INCREMENTAL  AND  DIRECT  DETECTION  PROBABILITIES 


Taking  the  assumptions  on  signal  and  noise  ampli- 
tudes stated  in  the  introduction,  and  further  assuming 
that  the  signal  only  need  exceed  the  noise  for  a detec- 
tion to  occur,  the  probability  of  detection  of  an  event 
phase  at  a station  is  given  by  a normal  distribution 
function : 


pD(m) 


(1) 


where  y is  the  mean  logarithm  of  noise  amplitude  at 
the  station,  a and  a are  the  standard  deviations  of 
the  logarithms  of  signal  and  noise  amplitudes  respec- 
tively, and  A is  the  calculated  amplitude  of  that 
phase  at  the  station  from  an  event  with  operational 
magnitude  m.  Computation  of  PD(m)  can  be  extended  to 
a network  if  desired  (Wirth,  1970).  Blandford  and 
Wirth  (1973)  have  analyzed  the  case  of  a station  or  a 
network  with  automatic  power  or  F detectors  and  have 
justified  the  use  of  (1)  for  those  stations.  The 
detection  curve  for  a given  station- region  pair  should 
take  the  form  of  equation  (1)  when  the  method  of  direct 
threshold  determination  is  employed.  Now,  when  a 
station  is  used  to  observe  events  for  a period  of  time, 

, A 

data  on  numbers  of  events  versus  observed  magnitude  m 
at  that  station  may  be  accumulated;  and  incremental 
thresholds  can  be  determined  from  this  data  by  maximum 
likelihood  or  least-squares  (Kelley  and  Lacoss,  1969). 


However,  these  thresholds  arc  biased  due  to  the  fact 
that  they  are  based  on  observed  magnitudes,  not  opera- 
tional magnitudes.  We  generally  need  to  know  P^(m), 
whereas  the  incremental  method  determines: 


PjCm) 


* 


logA  - 

o 

n 


u 


(2) 


A 

where  A is  now  the  estimated  amplitude  associated  with 

A 

m which  is  m plus  a random  error  variable  normally 
distributed  with  zero  mean  and  variance  due  to 
path-station  effects.  The  bias  of  the  incremental 
threshold  magnitude  relative  to  the  direct  one  can  be 
explicitly  calculated  as  follows.  First  consider  the 
number  of  events  detected  at  both  m and  m: 


Nj(m)  = N(m)*$ 


N^Cm)  = N(m)“l< 


m - y ' 


(o2  * a2)1'2 

v c n J 


m - * 


A constant  equal  to  the  distance  correction  factor  in 
magnitude  calculation  for  a particular  station-region 
pair  has  been  absorbed  into  y to  get  y'.  Taking  the 
fraction  of  events  detected  in  each  case  to  be  a percent 
and  equating  the  detection  probabilities  gives: 


m - y ' 
a 


(a2  * o2)1/2 
s nJ 


= 4* 


m 


a 


n 
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where  and  m^  are  the  operational  and  observed  magni- 
tudes at  which  a percent  detection  occurs.  Equating 
the  arguments  to  in  the  above  $ functions,  which 
are  equivalently  the  unit  normal  density  function 
integrated  to  2^,  and  eliminating  \i'  results  in: 


A 


m 


a 


* o2)V2 
n' 


(3) 


This  relation  gives  the  biased  observed  threshold 

magnitude  ma  which  is  determined  by  the  incremental 

method  when  the  operational  threshold  lies  at  m . 

a 

For  example,  if  the  operational  threshold  is  chosen 
to  be  at  90%  and  as=an=.3,  then  m g^m  g-0.16.  The 
bias  is  symmetric  about  the  501  threshold;  thus  for 
the  101  operational  threshold,  m ^m  1 + 0.16.  Table  I 
lists  che  10%,  501,  and  90%  threshold  magnitude  bias 
for  various  values  of  ag  and  on.  Note  that  the  above 
discussion  can  be  extended  to  true  and  observed  magni- 
tudes if  os2  is  replaced  by  o2  + a2  where  a2  is  the 
variance  of  the  source  effects  about  true  magnitude, 
Since  use  of  observed  magnitudes  gives  a biased  esti- 
mate of  the  operational  or  true  90%  threshold  but  an 
unchanged  estimate  of  the  50%  threshold,  in  practice 
better  comparisons  of  incremental  thresholds  among 
stations  might  well  be  made  at  the  50%  level. 

The  threshold  magnitude  bias  produced  by  a net- 
woi  of  stations  cannot  be  expressed  in  a simple 
equation  such  as  that  for  a single  station  (equation  3). 
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However,  if  we  can  estimate  numerically  what  the 
incremental  detection  curve  on  observed  magnitude  for 
a network  should  be,  then  by  comparison  with  the 
incremental  curve  on  operational  magnitude,  we  can 
determine  network  threshold  magnitude  bias.  Let  us 
begin  to  estimate  the  incremental  network  detection 
curve  with  an  expression  from  Herrin  and  Tucker  (1972) 
for  the  probability  density  ot  observed  m at  a station 
for  the  subset  of  events  with  operational  magnitude 
m detected  at  that  station: 


where  a replaces  (o^  + o^)1//2.  This  differs  from  the 

probability  density  of  all  m at  a station  for  all  the 

operational  magnitude  m events  (neglecting  the  fact 

that  m cannot  be  calculated  for  undetected  events) 

which  must  simply  be  normal  with  mean  m and  variance 

a2.  We  need  the  probability  density  of  m averaged  over 
s 

those  stations  which  detected  each  of  the  magnitude 
m events.  This  is  a sum  of  probability  densities 
similar  to  (4)  for  all  the  cases  of  "k  of  n"  detection, 
each  suitably  weighted  by  the  probability  that  exactly 
k of  n stations  detect. 

Unfortunately,  we  were  able  to  derive  the  distri- 
bution f (m)  only  for  the  case  of  two  stations  of  equal 
detection  capability;  this  derivation  is  given  in 
Appendix  1,  and  the  resulting  expression  (A12)  is: 
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Although  we  will  not  demonstrate  it  here,  the  expected 
value  of  m with  the  probability  density  given  by  (5) 
is  the  same  as  the  expected  value  of  m with  the  prob- 
ability density  given  by  (4),  in  other  words,  on  the 

A 

average,  the  expected  observed  magnitude  bias,  m-m, 
for  an  event  of  given  m will  be  the  same  regardless  of 
whether  one  or  several  stations  (of  equal  capability) 
detected  it.  (Of  course  the  bias  does  vary  as  a 
function  of  m.)  Herrin  and  Tucker  (1972)  have  demon- 
strated this  fact  without  having  to  derive  the  distri- 
bution of  f (m) . 

To  calculate  the  number  of  events  with  magnitude 
m that  will  appear  for  a two-station  network  or  a 
single  station,  we  integrate  f(m)  or  f(m),  respectively, 
weighted  by  the  number  of  events  detected  at  m,  over 
operational  magnitude  m: 

Nj(m)  = f (m) *PD(m) *N(m) .dm  (single-station) (6a) 
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Nd(m}  = /.co  [ f (m)  - Bj  >2  (m)  + f(m).B2^2(m)J  N(m)-dm  (6b) 

(two-station  network) 

wlicre  B^2  and  arc  the  binomial  probabilities  of 

one  and  two  station  detection,  respectively,  given  by: 

u _ r nA n n-k 
k , n ~ <h)P  (1 

where  p is  given  by  equation  (1)  and  q is  1-p.  i\c 
will  need  for  our  derivation  a relation  between  the 
number  of  events  which  occur  and  the  magnitude;  this 
is  commonly  expressed  as: 

log10N(m)  = -am  + 0 

or  as 


where  3 = ldp  and  a'=2.503a.  Kc  performed  the  inte- 
gration of  (6a)  and  ( 6b ) numerically  for  the  ease  of 

°n  * ^ * as'*5»  a=0.8,  0 =1.7x10^,  over  a finite 

range  of  m to  get  th  a network  incremental  curve 
and  a single-station  incremental  curve.  The  results 
arc  shown  in  figure  1 along  with  the  assumed  seismicity 
magnitude  relation  for  operational  magnitude  given 
by  1-0.8,  3 =1.7x10**.  In  this  figure  the  seismicity- 
magnitude  relation  is  plotted  as  a function  of  opera- 
tional magnitude  m,  while  the  incremental  curves  arc 
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plotted  against  observed  magnitude  m.  An  incremental 
threshold  determination  for  the  single-station  case  is 
made  by  observing  the  percentage  falloff  in  detection 
of  the  number  of  possible  events  predicted  by  fitting 
a straight  line  to  the  single-station  points  at  higher 
magnitudes.  The  single-station  50%  and  90%  threshold 
magnitude  on  operational  magnitude,  as  calculated  by 
equation  (1),  are  indicated  for  reference.  Note  that 
the  single-station  bias  for  the  90%  threshold  magnitude 
is  toward  a lesser  magnitude  as  predicted  by  equation 
(3)  and  that  the  bias  for  the  50%  threshold  magnitude 
is  zero  as  also  predicted  by  equation  (3).  The  hump 
in  the  network  incremental  curve  near  the  threrhold  is 
characteristic  and  is  confirmed  by  points  generated 
in  a simulation  of  the  two-station  network  by  a com- 
puter program  which  we  will  discuss  later.  This  hump 
is  an  unfortunate  happenstance  since  it  will  degrade 
our  capability  to  make  threshold  determinations  for 
networks  by  fitting  a straight  line  to  the  observed 
seismicity-magnitude  data,  Nj(m).  Consequently,  we 
cannot  easily  determine  the  threshold  bias  for  a net- 
work because  of  the  uncertainty  about  the  true  asymptotic 
slope  of  the  seismicity-magnitude  relation  when  given 
a d;splay  of  network  data  such  as  that  in  Figure  1. 

When  the  number  of  detected  events  Nj(m)  is  plotted 
against  observed  magnitude,  then  the  asymptotic  portion 
of  this  curve  where  Nj(m):sN(m)  is  displaced  from  the 
seismicity-magnitude  relation  plotted  against  true  or 
operational  magnitude  as  in  Figure  1.  This  displacement 
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is  due  to  the  fact  that  more  events  from  a lower  opera- 
tional magnitude  scatter  into  a given  observed  magnitude 
than  from  a higher  operational  magnitude,  be  assume 
that  the  observed  magnitudes  arc  equal  to  an  operational 
magnitude  plus  a random  variable  normally  distributed 
Wlth  zero  mean  and  variance  o‘.  We  can  calculate  the 
displacement  of  the  asymptotic  portion  of  the  observed 
seismicity-magnitude  line  relative  to  the  real  one  by 
finding  the  number  of  events  with  observed  m. 


/s  / oo 

N(m)  = f_( 


(2r) 


exp 


. N (m) *dm. 


Substituting  (7a)  in  the  above  gives 

Q A fCn 

N(m)  = YTl — /-<*>  C>P 
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which  is  easily  integrated  to  give: 


N (m) 


-a"m  a"^o^/2 
3 'c  c 


(7b) 


Let  m and  mR  be  the  observed  and  operational  magnitudes 

at  which  k events  appear;  then  setting  N(mk)  equal  to 

N(m.)  and  taking  the  natural  logarithm,  we  have: 

K ~ 


The  horizontal  displacement  of  the  asymptotic  portion 
of  the  seismicity-magnitude  line  is  just  a'vj  2,  which 
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f 


.should  he  roughly  .1  to  .2  magnitude  units  for  a single 

-station  for  reasonable  values  of  a and  aI 2.  We  see  in 

figure  1 that  the  single-station  observed  seismicity 

line  has  been  displaced  relative  to  the  real  one  .10 

outward  along  the  magnitude  scale.  The  displacement 

of  the  line  vertically  leads  to  an  overestimation  of 

seismicity  at  a given  magnitude,  given  by  the  factor 
1 n ‘ 31  0 s / ? i • , 

•/-,  which  is  approximately  1.15  for  the  single- 
station  case  in  Figure  1.  For  a network,  a2  is  divided 
by  the  number  of  stations  over  which  magnitude  is 
averaged  as  reflected  in  the  lesser  displacement  of 
the  two-station  line  in  Figure  1;  and  the  displacement 
of  the  observed  network  seismicity-magnitude  line  from 
the  true  one  should  be  insignificant  for  a network  of 
more  than  a few  stations. 


I he  above  derivation  on  displacement  of  the  seis- 
micity line  can  likewise  be  made  for  operational  versus 
true  magnitudes  simply  by  replacing  by  o^.  However, 
if  there  is  no  source  bias,  i.e.  am=0,  then^he  seis- 
micity line  based  on  operational  magnitude  will  not 
dilfer  from  that  based  on  true  magnitude,  and  the 
obsc-vcd  network  seismicity  line  will  approach  that 
based  on  true  magnitude. 


In  this  section  we  have  shown  that  thresholds 
(other  than  505)  determined  incrementally  for  a single 
station  will  differ  from  thresholds  based  on  opera- 
tional or  true  magnitude,  that  a hump  will  appear  in 
the  incremental  curve  for  a network,  thus  making 
threshold  determination  in  the  network  case  difficult, 
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and  that  seismicity  can  be  somewhat  overestimated  by 
incremental  detection  plots.  We  will  next  investigate 
some  of  the  problems  which  arise  in  connection  with 
the  use  or  cumulative  plots  to  determine  threshold. 
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CUMlihAT  1 VI;  AND  lNCRIiMLNTAL  Dl.Tl.CTION  PROBAIU L 111  IiS 

Often  it  seems  necessary  to  determine  thresholds 
by  plotting  the  cumulative  numbers  of  events  with 
magnitudes  greater  than  a given  magnitude  instead  of 
the  incremental  numbers  as  in  the  previous  section 
Clearly,  this  method  utilizes  the  same  data  as  the 
incremental  method  only  in  a modified  fashion.  As 
Lacoss  (1972)  has  already  pointed  out,  cumulatively 
determined  threshold  magnitudes  will  always  be  lowrer 
than  incrementally  determined  ones,  he  can  derive 
this  bias  analytically  as  follows,  first,  we  express 
the  incremental  and  cumulative  curves  as: 

l’jCm)  = N'd(in)/N'(m)  (9) 

i’cO")  = /-  Nd(m')dmV/-  NfiVjdm"  (10) 


If  we  assume  the  common  form  for  the  seismicity- 
magnitude  relation  given  by  equation  (7b),  then  P^(m) 
can  he  incorporated  in  Pj.(m) 


PcCm)  = j:  C'a'm'  I’,  ( m " ) dm'/ c'a'n'dZ'. 


\ou  assuming  P j (m)  to  be  the  normal  distribution  func- 
tion given  by  (2)  (strictly  speaking,  we  have  seen  that 
this  is  valid  only  for  stations,  not  networks)  and 
integrating  the  denominator  of  the  above  equation, 
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we  have: 
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A linear  trails  format  ion  of  the  variable  of  integration, 

A A 

m'  = m"  + p', results  in 
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he  now  make  the  substitution  of  the  irror  function 


n a a (m-p  ) -am  , aa 

l'c(m)  = -f  /j.^c  dm  * 


(11) 
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The  first  integral  in  (11)  is  well  known: 


I 


-a  m , A ^ _ 1- a (m-p  ) 
a c dm  = —re  v 

m-p  a 


The  second  integral  in  (11)  can  be  evaluated  (e.g., 
Abramowitz  and  Stegun,  1964,  p.  304)  and  reduces 
straightforwardly  to: 


— - — 


ihus,  substituting  these  in  (11)  and  performing  some 
algebraic  manipulation,  we  have: 
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Our  assumption  for  I5 j (m)  in  equation  (2)  requires 


PjCm)  • 4 


u n 

1 + erf 

’m  - y \ 

i^7/J 

Our  result  then  for  the  probability  bias  in  the  cumu- 
lativc  detection  curve  versus  the  incremental  is: 


l’c(m)  - P j Cm)  =4  exp 


2 ,2 

^ o a 

a'(m  - y")  + -H- — 


n - \}'  + u'o1,1 

1 - erf 

n 

i 2i/Zo  j 

' n ' 

(12) 


Lacoss  (1972)  has  previously  expressed  this  bias  in  a 
form  equivalent  to  (12).  We  list  in  Table  II  the 
magnitude  bias  of  the  cumulative  threshold  determina- 
tion relative  to  the  incremental  determination  pre- 
dicted by  (12)  at  the  90?o,  50%  and  10%  thresholds 
for  various  reasonable  values  of  o and  a;  and  the 
bias  is  significant  in  all  eases.  The  bias  given  by 
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(12)  applies  strictly  to  a single  station  since  Pj(m) 
for  a network  has  the  hump  discussed  in  the  previous 
section  and  cannot  he  expressed  as  a normal  distribu- 
tion function.  However,  we  can  again  expect  the  bias 
to  be  toward  lower  magnitude,  by  some  unpredictable 
amount . 

Note  that  the  bias  given  by  (12)  is  for  the  cumu- 
lative curve  relative  to  the  incremental  curve.  If 
the  90%  threshold  obtained  from  the  incremental  curve 
is  already  biased  low  as  discussed  in  the  previous 
section,  the  cumulative  threshold  determination  will 
be  even  further  biased  low.  For  example,  if  o -a  = .3 
and  a=l  and  the  90%  threshold  is  considered,  then  from 
Table  II  the  cumulative  threshold  will  be  .21  lower 
than  the  incremental  threshold  which  is  already  biased 
.16  lower  than  the  threshold  based  on  operational 
magnitude  in  this  case.  This  is  a total  bias  of  -.37 
magnitude  unit,  a significant  amount  indeed.  Since 
there  is  no  incremental  bias  for  the  50%  threshold, 
the  cumulative  bias  only  need  be  considered;  for  our 
example  it  amounts  to  -.34  from  Table  II,  again  a 
significant  amount 

We  have  now  derived  the  errors  to  be  expected 
when  determining  the  threshold  magnitude  of  a station 
using  magnitudes  calculated  at  that  station,  and  have 
indicated  what  errors  to  expect  when  networks  are  con- 
sidered rather  than  a single  station.  This  knowledge 
allows  us  to  remove  the  errors  and  make  a valid  com- 
parison between  empirical  thresholds  and  theoretical 
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ones  based  on  operational  magnitude.  What  corrections 
are  necessary  if  we  determine  empirically  station  or 
network  thresholds  relative  to  an  independent  magnitude 
scale  (e.g.,  surface-wave  threshold  on  mfa)  and  then 
wish  to  compare  them  to  theoretical  ones  based  on 

operational  magnitude?  This  is  the  topic  of  the  next 
section. 


CONVERTING  DETECTION  THRESHOLDS  FROM  M to  m. 

s b 

Figure  2 shows  a hypothetical  but  typical  group 
of  data  with  M * s from  a network  whose  threshold  we 
wish  to  determine  and  m^ ' s from  the  same  or  an  inde- 
pendent network.  By  plotting  the  percentage  of  events 
at  given  in^ ' s for  which  a surface  wave  was  detected, 
one  can  simply  determine  a probability  of  detection 
curve  for  surface  waves  based  on  event  mfa . The  pro- 
blem is  to  convert  a surface-wave  threshold  on  m^ 
which  was  determined  directly  using  the  given  m^ ' s to 
an  equivalent  threshold  on  M , or  vice-versa.  The 
simplest  approach  in  transforming  LR  thresholds  from 
the  Mg  magnitude  scale  to  the  m^  scale  is  to  equate 
detection  probabilities  for  corresponding  points  on 
the  Ms_mb  line  which  fits  a set  of  empirical  data  from 
the  network  or  station  under  consideration.  This 
approach,  however,  leads  to  prediction  of  LR  thresholds 
on  m^  which  are  too  low  when  compared  to  those  actually 
determined  on  m^ . The  problem  is  indeed  more  complex, 
and  one  must  take  into  account  the  distributional  pro- 
perties of  the  data,  such  as  M and  m,  scatter  about 

s b 

the  Mg - line  and  increasing  numbers  of  events  with 

decreasing  magnitude,  as  shown  in  Figure  2.  The 

correct  transformation  of  thresholds  from  the  M to  m, 

s b 

scales  and  vice-versa  will  now  be  presented. 

The  detection  probability  of  surface  waves  con- 
ditional upon  mb  can  be  expressed  in  terms  of  that 
conditional  upon  M by: 
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(13) 


POnb)  = r.J CMs|mb) -PCMS)  “dMs 

This  is  a useful  relation;  that  is,  we  can  convert 
predicted  thresholds  for  surface  waves  in  terms  of  Mg 
into  terms  of  mb  to  compare  with  empirical  thresholds 
determined  directly.  To  evaluate  (13)  in  practice, 
we  need  explicit  expressions  for  the  two  factors  in 
the  integrand.  The  second  factor  is  given  for  a 
single  station  by  equation  (1).  (For  a network,  (1) 
may  be  a satisfactory  approximation,  since  the  detec- 
tion probability  curve  for  a network  is  an  F distri- 
bution (Herrin  and  Tucker,  1972)  which  approaches  a 
normal  distribution  as  the  number  of  stations  in  the 
network  increases.)  The  first  factor  in  (13)  can  be 
expressed  by  the  basic  probability  identity: 

f(Ms|mb)  = f(mb,Ms)/f(i%)  (l4) 

The  joint  density  function  f(mb,Ms)  of  the  operational 
magnitudes  reflects  the  scatter  in  a typical  Mg  versus 
mb  plot;  this  density  function  is  expressed  by  intro- 
ducing the  true  magnitudes,  Mg  and  mb , each  pair  of 
which  are  linearly  related  in  a precise  manner  thus: 

M = am,  + b . ^ ^ ) 

s b 

We  assume  that  Mg  and  mb  are  normal  random  variables 
with  means  Mg  and  inb  and  standard  deviations  and 
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O , . Then  the  joint  density  function  of  the  complete 
mb 

population  of  events  arising  from  events  with  true 
magnitude  m^  is: 


i -Q(m.  ,M  ) .. 

<•(■%, m ) - r„  2,a  g-,-  ° ’ f(mb)-dmb-  (16) 

u ms  mb 


The  exponent  Q in  the  above  is 


<KvMs> 


(17) 


M could  just  as  well  serve  as  the  variable  of  integra- 
tion  in  (16)  . Now  let  the  form  of  the  density  function 
of  m^  be  given  by  equation  (7a)  . Then  it  is  apparent 
from  a comparison  of  equations  (.a)  and  (7b)  that  the 


density  function  for  mb  is: 


~ <*2  2 / 7 

a bmb  -a  bamb/2 
e 


Note  that  f(m^)  will  be  equivalent  to  f (mb)  only  if 
omb  = 0 or  a=0 . This  f(inb)  is  not  a proper  probability 
density  expression  in  that  it,  as  well  as  f(mb)  an^ 
thus  f(mb,Ms)  above,  does  not  integrate  to  unity; 
however,  the  result  we  need,  f(Mslmb)»  is  a ratio  of 
these  as  given  by  (14)  and  is  finite  and  intcgrable 
to  unity.  Using  (15)  in  (17)  and  then  substituting 
(17)  and  (18)  in  (16)  and  performing  the  integration 
results  in 
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ctV'V  = 


exp 


6 b °xp  — 2 


A 2 
a b°mb 


[Mg  - (amb  + b)] 
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~1  7 2 2 

°ms  + 3 °mb 


*1 
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(19) 


Using  (7a)  and  (19)  in  (14)  we  arrive  at 


f(Mjm.)  = 


rn  (o 


ms 


2 2.1/2 
a “mb5 


exp 


2(o2  + aV,) 

v ms  mb' 


2 2 

[Ms  - (a.b  * b “ahaomb) 1 

Note  that  this  is  a normal  probability  density  function 

with  a mean  a.  aa2,  units  below  the  true  M -m,  line  and 
_h  mb9  ? , , 

a variance  o + a o , . Equation  (20)  should  accurately 
ins  m b 

describe  the  density  of  points  along  any  vertical  line 

in  an  M -m,  plot  such  as  Figure  2,  provided  that  the 
s b 

1,R  phase  from  all  events  with  the  given  mb  is  detected 
by  the  station  or  network.  To  illustrate  (20)  we  take 
a simple,  but  realistic,  case:  a=l  and  b=0  (Ms=mb) , 

N “ e"a  bmb  (natural  seismicity),  and  oms=<Jmb 
(equivalent  M and  mb  scatter).  The  conditional 
distribution  of  Mg  is  then: 


I (20) 
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f(>g,Jb)  = ■^777177—  c-'l>  (-7 


ms 


/ 

2' 

,4 

~1TTT, 

ms 

• _ 

J 

(21) 


So,  with  Ms«mb,  wc  have  a normal  density  function  for 
Ms  with  mean  and  variance  20;  ; the  mean  of 

the  distribution  lies  a'o^  below  the  true  M -m  line. 
Taking  the  above  case  again,  but  now  allowing  o ,=0 

(no  scatter  in  mb)  and  u'-O  (uniform  seismicity)1,  we 
get: 


f(M  | m.  ) = c\n 

1 

K ■ lnbl 

2- 

ms 

"2 

'ms  J 

Here  wc  have  simply  a normal  density  function  for  M 
with  mean  inb  and  variance  o“s.  This  is  the  density* 
function  proposed  by  Harley  and  llciting  ( 1972)  and 
Lacoss  (1971);  however,  (21)  is  the  proper  expression, 
and  it  includes  a downward  displacement  of  the  mean 
due  to  increasing  seismicity  with  decreasing  magnitude. 

he  now  know  both  factors  needed  in  (13)  to  per- 
form the  conversion  of  surface-wave  thresholds  from 

Ms  to  rnb  scalcs-  However , we  note  that  in  dealing 
with  observed  magnitudes  wc  need  P(mb)  rather  than 
P(mb).  If  wc  assume  that  mb  scatters  normally  about 
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mb  variance  osb,  the  form  of  the  above  expression 

is  preserved  with  (omb  ♦ replacing  omb  every- 

where. The  result  of  using  mb  rather  than  is  to 
predict  thresholds  by  (13)  which  lie  at  a lower  magni- 
tude than  observed;  that  is,  any  probability  calculated 
by  (13)  for  mb  really  applies  to  a higher  mb . If  the 
mb  is  a network  average,  the  overall  effect  is  small, 
less  than  .1  magnitude  unit  in  the  threshold  value. 

A similar  consideration  of  N1  lather  than  M would 

s s 

provide  the  necessary  distribution  f(Ms,mb)  if  one 
desired  to  use  Mg  threshold  calculations  given  by  (2). 

It  is  now  possible,  using  (20)  and  (1),  to  evalu- 
ate (13),  and  it  is  seen  that  this  evaluation  requires 
knowledge  or  estimation  of  several  parameters.  In  our 
application  of  (13),  we  therefore  estimate  an  empirical 
Ms-mb  correct  the  Ms"m^  line  for  the  Herrin  and 

Tucker  bias  at  low  magnitudes  (as  discussed  in  a later 
section),  choose  a seismicity-magnitude  relation  from 
some  independent  source,  and  assume  a normal  distribu- 
tion for  the  magnitude  scatter,  the  variance  of  which 
depends  on  source  bias  and  path  effects. 

A biased  estimate  of  M by  the  station  or  network 

due  to  favorable  or  unfavorable  travel  paths  for  LR 

affects  the  surface-wave  detection  threshold  on  the  mL 

b 

scale,  calculated  by  (13),  through  the  empirical  Ms~mb 

relation.  Empirical  surface-wave  thresholds  on  mL  will 

b 

be  equally  affected  though,  and  so  this  $s  bias  need 
not  be  considered  when  comparing  calculated  with 
empirical  detection  thresholds  on  mb#  Similarly,  any 
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bias  in  values  can  be  ignored  in  converting  surface- 
wave  thresholds  from  the  M.  to  mb  scales,  because  the 
empirical  scale  again  compensates  for  this  bias. 

However , in  reference  to  other  stations  or  networks, 
the  m^  bias  cannot  be  ignored.  The  reason  is  that, 
if  one  were  to  make  comparisons  of  LR  thresholds  on 

A _ 

m,  between  two  stations  or  networks  whose  sources  of 
*b  , . 

m,  values  were  not  identical,  then  discrepancies  could 
b 

arise  if  one  set  of  n^'s  were  biased  relative  to  the 
other.  These  discrepancies  can  only  be  eliminated  by 
estimating  and  removing  the  relative  bias  in  observed 

mb* 

The  calculated  detection  curves  in  Figure  3 show 
that  the  effect  in  (13)  of  using  a=l,0  (equation  21) 
is  to  raise  the  calculated  90%  threshold  on  the  mb 
scale  by  an  additional  .1  to  .2  magnitude  unit  over 
the  ease  where  a=0.0  (equation  22)  when  a normal  dis- 
tribution with  ams  = G,mk=  • ^ assume(^  f°r 

remainder  of  the  detection  curve  is  also  shifted  to 
higher  mb,  though  not  by  equal  amounts.  This  is  in 
agreement  with  simulation  results  to  be  given  later. 

One  may  want  to  perform  the  inverse  of  the  above 
problem;  that  is,  one  may  have  determined  empirically 
surface  wave  detection  probabilities  on  the  mb  scale 
and  desire  to  transform  them  to  the  Mg  scale.  The 
apparent  (but  incorrect)  procedure  is  simply  to  sub- 
stitute M for  m,  and  vice  versa  in  equation  (13)  and 
s b 
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proceed  with  the  numerical  integration  as  before,  only 
now  to  obtain  detection  probabilities  at  the  incremental 
Ms  s.  This  procedure  is  invalid  since  the  empirical 
probability  of  detecting  a surface  wave  from  an  event 
of  magnitude  mb  is  a result  of  looking  for  the  surface 
waves  from  events  of  varying  Mg  at  that  mb  , not  just 
a particular  Mg . Therefore  one  must  assume  some  detec- 
tion probability  curve  as  a function  of  M in  order  to 
evaluate  P(mb)  which  would  be  within  the  integral  when 
Ms  and  mb  were  interchanged;  but  this  detection  prob- 
ability curve  P(Mg)  is  exactly  what  we  are  attempting 
to  determine. 

We  conclude  that  there  is  no  direct  way  of  going 
from  an  mb  to  a Mg  scale  for  surface  wave  detection 
probabilities.  This  suggests  than  an  iterative  algo- 
rithm should  be  used  to  arrive  at  the  correct  detection 
probability  curve  on  the  Mg  scale.  To  do  this  we  take 
a starting  P(Mg)  curve  and  perturb  its  mean  (501  point) 
and  variance  while  making  the  transformation  according 
to  equation  (13).  This  requires  that  an  Ms-mb  relation- 
ship for  the  network  or  station  under  consideration, 
a seismicity-magnitude  relationship,  and  a standard 
deviation  for  mb  and  Mg  error  be  chosen  in  order  to 
calculate  the  factors  in  (13).  The  calculated  curve 
on  the  mb  scale  P(mb)  which  best  approximates  the 
empirical  one  indicates  the  correct  M detection  prob- 
ability curve,  P(Mg).  Figure  4 gives  the  results  of  a 
program  which  performs  this  algorithm,  in  this  case  for 

V1*0  mb * a=1*°»  and  amb=ams  = 0* 25  for  b°dy-wave  magnitude 
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and  surface  wave  magnitudes,  probably  a conservative 

value.  In  the  same  figure  the  curve  resulting  when 

M =1.2  m.-l.O  is  used  demonstrates  the  sensitivity  of 
s b 

t!ie  spread  of  the  calculated  I’(MS)  curve  to  the  slope 

of  the  M -m,  relation.  The  resulting  P(M  ) curve  is 

s b s 

also  dependent  on  the  intercept  of  the  Ms - relation 
and  would  slide  along  the  horizontal  axis  according 
to  the  exact  value  of  this  intercept.  We  also  show  in 
Figure  4 the  result  when  a=0.0  is  used;  evidently  the 
threshold  magnitudes  on  would  be  significantly  over- 
estimated if  seismicity  constant  with  magnitude  were 
assumed  when  the  true  relation  had  a=1.0. 

This  section  has  provided  a procedure  for  trans- 
forming thresholds  from  one  magnitude  scale  to  another, 
using  LR  thresholds  on  m^  and  as  the  example.  Again, 
the  results  are  strictly  for  a single  station; 
results  for  some  networks  may  be  closely  approximated 
by  those  for  a single-station  though.  The  importance 
of  the  seismicity-magnitude  relation  in  the  threshold 
theory  was  again  demonstrated;  and  in  this  case  of 
transforming  thresholds  between  magnitude  scales,  the 
importance  of  the  functional  relation  between  magni- 
tudes on  the  two  scales  became  evident. 
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DETERMINING  THE  TRUE  Mg-mb  LINE 
FROM  OBSERVED  DATA 


One  of  the  important  problems  in  seismic  discrimi- 
nation is  to  estimate  the  Mg -mb  relation  for  a given 
group  of  observed  pairs.  Determination  of 

the  true  relation  at  or  below  the  threshold  for  sur- 
face-waves (usually  found  to  be  higher  than  the  body- 
wave  threshold)  is  complicated  by  an  upward  Mg  bias 
in  the  observed  data  points  for  this  range  as  discussed 
in  Herrin  and  Tucker  (1972)  and  earlier  in  this  report. 
It  is  possible,  however,  in  theory  to  reconstruct  the 
true  M -mb  relation  from  the  observed  data  although 
there  will  be  difficulties  in  practice.  The  procedure 
is  as  follows: 

For  simplicity,  let  Ms«mb,  oms=omb,  and  ass=osb- 
Results  for  the  general  case  where  parameters  vary 
can  be  apparent  from  the  results  we  will  attain  with 
these  simplifications.  Equation  (21)  gave  the^density 
function  f(Ms|mb)  for  a single  station;  for  f(Ms|mb) 
appropriate  to  observed  magnitudes,  we  add  the  effects 
of  signal  variance  to  get: 


f(Ms lmbJ 


/ 


exp 
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\ 
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The  mean  of  this  normal  density  function  parallels  the 
true  M -m,  line,  and  t li e variance  is  constant  over  all 
mb.  Then,  for  a given  e,  0 < e < 1,  the  integral  of 
the  density  function  at  a given  in^  : 


"b . > d”s 

1 


(23) 


must  be  satisfied  with  an  M such  that  all  (M  ,m,  ) 

b . v s . * b .* 

pairs  will  parallel  the  trueiM.-m,  line.  One  cJn  1 

S D /v 

specify  exactly  that  z which  will  cause  Mg  to  fall 
on  the  true  M -m.  line.  In  the  case  M =m,  J let  the 
integral  (23)  be  evaluated  to  =mb 


M 


m,  - ct.'  (a c + o “ ) 
b^  b v ms  ssJ 


-i  \ 
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b. 


(2)  1/2(o’  + a2  )l/Z 

K J v ms  ss; 


dM  * 
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J 

Using  the  properties  of  the  normal  density  function, 
this  can  be  transformed  to  an  integral  over  the  unit 
normal  density: 


= / 


2 x 2 ,1/2 
a,  (a  +o  ) 
b v ms  ss 


(2) 


T7T 


Z( z) dz. 


(This  result  would  be  unchanged  if  Ms=amb+b.)  For 
example,  if  = 1 and  °ms=0ss= • 21 2 , then  z Z .32 
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r--. 


would  be  required  to  have  M on  the  M -m.  line: 

s • s b * 

repetition  of  this  at  varioui  m,  would  define  the 

b . 

true  line  by  a series  of  piints.  This  procedure 

could  be  applied  only  in  the  ideal  case  where  both 
surface  and  body  waves  were  detected  at  the  statir 
for  all  events. 

With  real  data,  one  must  account  for  thresholds; 
that  is,  all  the  points  predicted  by  the  above  density 
function  will  not  appear  in  the  data  due  to  nondetec- 
tions. The  observed  density  of  points  is  related  to 

f(Mjmb  ) by 
i 


f<Ssl”b.>  = W'l  V^^sw^V  • pbw<”b.)] 


where  Psw  and  P^  are  the  detection  probabilities  for 
surface  and  body  waves  in  the  form  given  by  equation 

A A 

(2)  earlier.  By  definition  here,  f (M°  | m^0)_<l . Using 
this  reduced  density  function  for  observed  data  and 
substituting  in  (23),  we  have  then: 


(24) 


A /\ 

Using  an  actual  plot  of  observed  (M  ,m,)  pairs,  one 

can  evaluate  (24)  numerically  down  to  the  M at 

s • 

which  the  integral  equals  e.  In  order  to  doHhis,  the 

A 

seismicity-magnitude  relation  log  N = -a,  m,  + B and 

b b 
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P and  P,  must  be  estimated  by  standard  procedures, 
sw  bw 

From  the  seismicity-magnitude  relation,  the  number  ot 
events  expected  to  lie  at  each  mb  , say  NL  , can  be 
determined ; and  if  we  let  N.  be  tAe  number  of  observed 

J /V 

events  in  each  increment  of  M at  that  m,  , we  can 

s u j_ 

approximate  equation  (24)  by: 


or,  removing  the  constants  from  the  summation, 

A 

M 


NiPbwCmb  )c  “ - l: 


M.  =M _ 
bj  bmax 


N /P  , ,(MC  ) 
y swv  Sj 


(25) 


where  M.  is  the  largest  M aeon  for  the  given  mb  . 

smax  ' i 

A ^hoice  of  c can  still  be  made  as  shown  above  to  make 

computed  M 's  lie  on  tlie  true  Ms-mb  line.  But 
in  most  cases  it  would  be  preferable  to  use  a smaller 
c,  say  0.1  or  0.2,  so  that  the  method  can  be  applied 
down  to  a lower  mb ; for  with  e larger,  one  simply 
runs  out  of  events  due  to  nondetection  of  surface 

A 

waves  before  the  summation  approaches  NiPbw(,nb  . ) e • 

An  example  of  the  application  of  equation  (25) 
will  be  shown  in  the  section  on  simulation.  The 
procedure  suggested  here  is  valid  only  for  single- 
station M -mb  plots.  It  requires  ather  precise 
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estimation  of  the  probability  detection  curves  of  the 
station  for  body  and  surface  waves;  and  if  such 
precision  is  unattainable,  then  one  should  not  extend 
the  procedure  below  roughly  the  expected  90%  threshold 
for  either  surface  waves  or  body  waves. 

Before  proceeding  to  a simulation  approach  to 
examine  network  threshold  problems,  we  will  tieat  in 
the  next  section  one  additional  effect  on  threshold 
whose  simple  aspects  can  be  handled  analytically  - the 
matter  of  correlated  signals  and  noise  among  stations. 
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CORRELATED  SIGNALS  AND  NOISE 

Noise  background  at  several  stations  of  a 
network  cannot  be  assumed  to  fluctuate  independently 
because  of  seismic  events  or  swarms  of  events  which, 
considered  as  noise,  can  produce  temporarily  higher 
background  at  two  or  more  stations  and  because  of 
mi c rose ismic  origins  which  can  simultaneously  affect 
many  stations  within  range  of  these  disturbances.  In 
applying  equation  (1)  to  a network  (Wirth,  1970)  we 
must  account  for  correlation,  if  any,  among  the  noise 
amplitudes  for  the  stations  in  the  network.  The 
effect  of  correlation  in  general  depends  on  the  magni- 
tude distribution  of  the  population  of  events  under 
consideration „ 

If  due  to  source  bias  a set  of  events  have 
operational  magnitudes  varying  about  their  true  magni- 
tude, then  for  this  set  of  events  there  will  be  a 
correlation  between  the  magnitudes  observed  at  dif- 
ferent stations.  In  other  words,  source  bias  implies 
signal  correlation  in  reference  to  true  magnitude. 

If  the  noise  is  also  correlated  as  described  above, 
it  is  of  interest  to  determine  the  correlation  be- 
tween stations  of  the  arguments  of  expression  (1). 

The  mathematical  expressions  of  these  relationships 
will  be  useful  in  assessing  the  effects  of  correlation 
on  thresholds.  We  define  the  observed  magnitudes  m. 
and  m2  at  two  stations  for  events  of  true  magnitude 
m as  the  sum  of  two  normal  processes: 
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m1  = G(m,om)  + G1(0,os) 
m2  = G(m,om)  + G2(0,os) 


Here  G is  the  random  normal  variate  representing  scatter 
of  operational  magnitude  about  the  true  magnitude,  and 
Gj  and  G2  are  random  normal  variates  representing  path 
effects  which,  when  added  to  the  operational  magnitude, 
produce  observed  magnitudes.  The  correlation  p~ 
between  and  m^  is  given  by: 


m 


o 

m 

“T~  1 

o +o 
m s 


(26) 


Note  that  p~=0.5  if  a =a  /O  and  p~  = 0 if  a =0.  If  m- 
m s m mm  1 

and  m?  are  normal  processes  with  correlation  as  in 

z 22  1/2 
(26),  mean  m,  and  standard  deviation  (om  + og)  and 

if  and  u'2  arc  zero-mean  normal  processes  with 

parameters  p and  a ; the  correlation  p between 

(™l"u>l)  ant*  ^m2"y>2^  8^ven  b/: 


P = 


7 ? 2 

p~(o^  + a ) + p o 
^mv  m s'  Kn  n 

r~  t~  2 

a + a + a 
m s n 


(27) 


Note  that  p=p*=p  if  p*  and  p are  equal, 
m n m n 
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Ivc  can  obtain  analytical  predictions  for  the  effect 
of  correlation  on  a network  threshold  in  only  the 
simplest  of  cases.  A two-station  network  with  both 
stations  having  an  equal  probability  of  detection  of 
events  can  be  treated  through  the  bivariate  normal 
distribution.  Consider  the  case  of  random  normal  noise 
and  signals  so  that  the  correlation  of  the  variables 
ni  1 “ P 1 *'md  in  2 - M ^ 2 *s  as  «ivcn  equation  (27). 

Bivariate  normal  curves  for  p=  0.0,  0.5  and  1.0  were 
used  to  obtain  the  ''one  or  more"  detection  probabilities 
for  this  hypothetical  two-station  network  based  on 
p =5.00,  (o2  ♦ op^2  = 0.0  and  an=0.2  as  sho^n  in  Figure  5. 
The  effect  of  perfectly  correlated  noise  is  to  raise 
the  50°o  and  90°o  thresholds  by  little  more  than  0.1 
magnitude  unit. 

Ivc  cannot  predict  analytically  tl  e noise  correla- 
tion effect  for  realistic  networks,  though,  and  a 
simulation  experiment  is  required,,  A program  was  written 
to  generate  populations  of  signals  and  noise  at  the 
various  stations  of  a network  according  to  the  desired 
statistical  parameters  and  to  select  out  those  "events" 
which  were  detected.  Tl'  .s  program  will  be  described 
in  more  detail  in  a later  section.  To  check  the  pro- 
gram, we  compared  its  detection  curve  in  the  "one  or 
more"  detection  case  for  a two-station  network  with 

P '=  5 . 00  , o =0.2,  and  p=.95  for  the  noise  and  o2=o2=0 
* * ms 

for  the  signals  to  the  theoretical  prediction  as  shown 
in  Figure  5.  (The  correlation  coefficient  in  the 
simulation  program  was  chosen  to  be  .95  since  exactly 
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1.0  would  give  negative  eigenvalues  for  the  noise 
matrix.)  The  simulated  and  predicted  curves  seem  to 
agree  well,  keeping  in  mind  this  slight  difference  in 
the  value  of  the  noise  correlation.  We  proceeded  next 
to  simulate  a "two  or  more"  detection  curve  for  the 
nine-station  Long  Period  Experimental  network.  The 
noise  amplitudes  and  noise  correlation  coefficients 
were  taken  from  von  Seggern  0974);  the  p averaged 
roughly  0.4„  An  epicenter  at  50N  and  120E  was  arbi- 
trarily chosen,  and  the  simulation  was  run  for  o =.2 
and  os=.3,  °m=0»0»  and  p*=0.0.  The  resulting  detection 
curve  as  a function  of  true  magnitude  is  shown  in 
Figure  6 along  with  the  one  theoretically  predicted 
for  p=0.  Correlated  noise  has  apparently  little  or  no 
effect  in  this  case. 

All  the  foregoing  discussion  on  effects  of  corre- 
lated noise  expresses  thresholds  in  terms  of  true 
magnitude  and  assumes  that  the  number  of  events  versus 
magnitude  is  uniform.  The  true  seismicity-magnitude 
relation  is  commonly  approximated  by  log1QN^amb+g 
though,  and  we  must  account  for  this  in  our  simulation 
experiment.  We  will  include  such  a case  with  correlated 
noise  on  a ten-station  network  in  the  next  section, 
which  will  describe  our  simulation  procedure  and  results. 
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SIMULATION  LXPLRI  MLNT 


NccJ  Cor  a Simulation  Approach 

While  we  have  been  able  in  the  previous  sections 
to  derive  analytically  many  useful  relationships  between 
the  thresholds  a:u  properties  of  a single  station,  or 
in  some  cases,  a two-station  network,  many  results 
required  for  practical  evaluation  of  a real  network 
could  not  he  obtained.  Lor  example: 

1.  A network  does  not  have  a response  similar  to  a 
single  station,  and  analysis  of  the  network  response 
becomes  prohibitively  complicated  when  the  different 
stations  of  the  network  have  noise  levels  with  differing 
means  and  variances  and  arc  at  different  distances  from 
the  epicentral  area. 

2.  In  a real  network  the  I1  threshold  influences 
the  LR  threshold  by  virtue  of  the  fact  that  one  does 
not  attempt  to  detect  LR  unless  I’  has  been  detected. 
While  this  will  not  raise  the  LR  threshold  if  the  P 
threshold  is  far  below  the  LR  threshold,  many  real 
networks  do  not  have  this  property.  lor  example, 
detection  of  LR  by  the  LPI  network  in  response  to  NOS 
epicenters  does  not  (von  Scggern,  1974) , and  a network 
of  LASA  stations  would  not.  We  shall  show  for  this 
latter  case  that  the  thresholds  arc  affected  to  the 
degree  of  0.1-0. 2 magnitude  units. 

3.  Source  bias  due,  for  example,  to  differing 
fault  mechanisms  of  an  earthquake  or  differing  depths 
of  burial  of  an  explosion  influences  thresholds. 
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4.  Correlation  of  mean-square  noise  amplitudes 
among  stations  influences  thresholds;  but  as  we  shall 
see  in  the  simulations  the  effect  is  almost  negligible 
for  realistic  correlation  levels  when  source  bias  also 
exists . 

5.  Incremental  and  cumulative  histograms  of  the 
number  of  events  as  a function  of  mb  for  which  LR  has 
been  detected  are  often  of  interest.  However,  analysis 
or  such  curves  has  thus  far  eluded  analytical  investi- 
gation even  for  single  stations. 

6.  Determination  of  the  true  Ms_mb  relationship 
by  fitting  straight  lines  to  observed  data  sets  is 
unreliable  due  to  magnitude  bias  introduced  around  the 
station  and  network  thresholds.  This  bias  can  be 
evaluated,  and  thus  corrected,  by  means  of  simulation. 

M.ThQRTll  Predictions 

As  we  shall  see  below,  it  has  been  possible  to 
use  the  standard  program  NhTWORTIi  (Wirth,  1970)  to 
compute  thresholds  for  true  magnitudes  if  there  is  no 
source  bias  or  noise  correlation,  and  to  compute 
thresholds  with  respect  to  operational  magnitudes  if 
there  is  no  noise  correlation.  Some  modifications  to 
Nil  TWO  RTH  are  required  in  order  to  compute  thresholds 
when  noise  correlation  is  present  or  to  compute  thresh- 
olds with  respect  to  true  magnitudes  when  signal  cor- 
relation (source  bias)  is  present.  The  modified  pro- 
gram, NLTWCORR,  does  not  compute  thresholds  with  respect 
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to  Observed  magnitudes.  Only  the  full  simulation 

program  MSBM.T,  to  he  discussed  below,  will  accomplish 
that . 

A detailed  discussion  of  the  new  XI.TKCORR  techni- 
ques is  contained  in  Appendix  II;  however,  a short 
outline  will  he  given  here.  The  user  of  NI.TIVCORR  now 
lias  the  option  of  reading  in  a signal  correlation 
matrix  and  a noise  correlation  matrix.  Iquation  (26) 
may  he  used  to  calculate  the  elements  of  the  signal 
correlation  matrix  from  a^given  variance  o2  of  source 
bias  and  given  variance  of  signal  amplitudes  for  an 
event.  At  each  magnitude  the  technique  outlined  by 
Shumway  and  Blandford  (1970)  is  then  used  to  generate 
sets  of  random  signal  and  noise  amplitudes  with  the 
given  correlation  structures.  hach  set  of  such  ampli- 
tudes represents  an  event  and  supplies  arguments, 
whose  correlation  coefficient  is  given  by  equation  (27), 
in  equation  (1)  for  all  stations  of  the  network.  If 
this  argument  is  greater  than  zero  a detection  by  that 
station  is  declared.  If  k or  more  stations  detect 
the  event,  then  this  increases  the  probability  of  a 
"k  or  more"  network  detection  for  an  event  of  the 
magnitude  under  consideration.  A sufficient  number  of 
random  events  is  generated  at  each  magnitude  to  produce 
a stable  probability  estimate. 

Table  III  describes  the  only  network  considered 
m detail  in  this  report.  it  has  been  designed  to  fit 
as  closely  as  possible  ten  l \SA-cquivalent  station3; 
each  of  v.iich  is  65°  from  the  Kamchatka- Kuri  1 region. 
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Ihe  y'  and  u'  v.lues  are  the  estimated  SOI  threshold 
magnitudes  for  LASA,  as  will  be  discussed  in  a later 
section.  The  estimates  for  the  signal  variances  are 
consistent  with  the  work  of  von  Seggern  (1972)  and 
1 ricsson  (19,1a),  and  is  consistent  with  short- 
period  values  for  LASA  beams  obtained  by  Chang  (1973). 
The  same  was  assumed  for  long-period  beams.  The 
a value  of  0.8  is  roughly  correct  for  the  Kamchatka- 
Kuril  region  earthquakes  on  the  mfa  scale  (W.  Dean  and 
R.  Abner,  personal  communication).  We  have  assumed 
Ms'1.0  mb*  ^ more  reasonable  formula  is  Ms=1.0  m^-0.8, 
calculated  from  a LASA-Kuriles  N^-m^  diagram  in  Mack 
(1971).  Then  if  an  observed  Mg  threshold  is  3.00,  one 
would  translate  that  into  3.80  for  comparison  with 
this  simulation  experiment.  The  reverse  procedure 
w o u 1 d need  to  be  carried  out  to  translate  simulation 
Ms  thresholds  downward  for  comparison  with  observazion. 

In  this  report  we  consider  only  networks  of 
type  I in  which  a station  does  not  attempt  to  detect 
LR  unless  it  has  detected  P.  A network  of  type  II 
uould  be  one  in  which  all  stations  look  for  LR  if  a 
specified  number  of  them  detect  P.  The  program  MSBNET 
may  be  run  for  either  type  of  network.  For  the  examples 
considered  in  this  report  the  changes  in  90$  thresholds 
are  less  than  0.1  between  runs  for  the  two  types. 

In  Table  IV  we  have  the  NETWORTH  and  NETWCORR 
results.  Figure  7 gives  complete  curves  for  the 
standard  NETWORTH  results.  The  90$  points  listed  in 
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Table  IV  show  that  there  is  approximately  0.2  magnitude 
difference  in  the  thresholds  between  the  standard 
NET WORTH  i?sults  and  the  NETWCORR  results  with  source 
bias  (with  and  without  noise  correlation)  if  the  true 
magnitudes  are  used  as  reference.  If,  however,  the 
operational  magitudes  are  used  as  reference  the  dif- 
ference is  closer  to  0.1  magnitude  units.  In  both 
cases  the  predicted  thresholds  with  source  bias  are 
higher  than  without.  The  differences  in  thresholds 
are  generally  smaller  but  still  have  the  same  sign 
for  the  503  and  103  levels.  Note  that  503  true  and 
operational  thresholds  for  the  case  of  source  bias 
present  are  nearly  the  same,  with  or  without  noise 
correlation. 

As  far  as  we  know  these  NETWORTH  results  are 
correct  if  the  underlying  reasonable  assumptions  are 
correct.  This  does  not,  however,  mean  that  they  will 
agree  with  observation  because  observed  magnitudes 
are  different  from  both  operational  and  true  magnitudes. 
Thus  NETWORTH  may  perfectly  predict  the  probability  of 
detecting  an  explo: ion  with  a specified  yield  (trans- 
lated into  true  magnitude)  and  yet  not  be  verified  by 
observation  if  observed  magnitudes  (instead  of  yields) 
are  used  to  determine  the  threshold.  The  only  way  in 
practice  to  establish  the  verifying  relationships  is 
by  means  of  the  simulation  techniques  to  be  discussed 
in  the  next  section. 
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MSBNLT  Simulation  Techniques 


A complete  discussion  of  program  MSBNET  is  given 
in  Appendix  III;  here  we  give  only  a short  outline,, 

The  user  specifies  station  parameters,  means  and 
standard  deviations  of  the  noise  at  the  stations,  the 
standard  deviation  of  the  signal,  the  seismicity  curve 
together  with  the  magnitude  limit  below  which  true 
magnitude  events  are  not  to  be  randomly  generated,  the 
portion  of  the  signal  variance  due  to  source  bias,  the 
noise  correlation  matrix,  and  the  magnitude  limits 
between  which  M -m^  regressions  are  to  be  calculated. 

Let  us  first  consider  the  case  where  there  is  no 
interstation  noise  correlation.  Consider  the  lowest 
0.1  true  magnitude  increment  in  which  random  events 

3 

are  generated.  For  a typical  run  it  will  comprise  10 
"events"  with  the  appropriate  true  m^  value  and  the 
corresponding  true  M.  value . Consider  now  the  first 
network  station.  To  each  of  the  k pail  of  true 
magnitudes  are  added  two  random  normal  numbers  with 
zero  mean  and  the  appropriate  standard  deviations.  One 
of  these  numbers  is  the  source  bias  and  for  each 
individual  event  will  be  the  same  at  all  stations. 

The  other  number  is  the  "unexplained"  signal  variance 
due  to  travel  path  effects,  and  this  will  be  different 
for  each  station  and  each  event.  Signal  amplitude 
correlation  in  this  situation  may  be  calculated  by 
equation  (26).  The  argument  of  equation  (2),  and  thus 
the  probability  of  detection,  may  then  be  calculated 
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for  both  P and  LR  using  the  newly  generated  observed 

A A 

and  Mg  values.  Two  uniform  random  numbers  are  then 
generated  in  the  interval  (0,1);  and  if  these  numbers 
are  less  than  the  corresponding  probabilities  of 
detection,  the  phase  is  "detected"  at  that  station. 

To  make  the  following  discussion  clear  it  seems 
necessary  to  be  somewhat  abstract.  For  ease  of  reference 
we  choose  to  use  the  same  variable  names  as  are  used 
in  program  MSBNET.  The  variable  names  are  suggestive 
of  their  role  in  the  program. 

Suppose  the  P wave  from  event  K is  detected  by  a 
station;  then  NOMBD(K)  is  increased  by  1,  and  the 
observed  value  of  mb  is  added  to  AMBM(K) . The  same 
procedure  is  followed  for  LR  and  Mg , using  NOMSD(K) 
and  AMSM(K)  if  the  P wave  is  detected  for  network 
type  I,  or  in  either  case  for  type  II.  After  all 
stations  for  all  events  for  all  magnitude  increments 
have  been  considered,  if  NOMBD (K) >NRMBD  then 
AMBN (K) /NOMBD(K) =mb  is  computed.  (NRMBD  is  the  number 
of  P-wave  detections  required  for  a network  detection; 
if  NOMBD(K) <NRMBD  we  go  to  event  K+l.)  Then  the 
counter  in  the  NMBNTD  appropriate  to  m^  is  increased. 

Then  if  NOMSD(K) >NRMSD,  we  compute  AMSM(K) /NOMSD(K) =Ag 
and  increase  the  counter  in  the  appropriate  NMSNTD 
magnitude  cell  and  in  the  NMBMSD  magnitude  cell  deter- 
mined by  m^.  (But  if  NOMSD(K) <NRMSD,  the  number  of  LR 
detections  required  for  a network  detection,  we  go  to 
event  K+l.)  After  all  events  are  considered,  regression 
curves  may  be  computed  on  the  surviving  m^ , pairs. 
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The  entire  procedure  may  then  be  repeated  for 
statistical  stability,  starting  with  a new  random 
number  and  accumulating  more  events  in  the  NMBNTD, 
NMSNTD,  and  N MB  MS  I)  arrays.  Similar  statistics  are 
kept  for  the  individual  stations.  Finally,  the  numbers 
required  for  the  incremental,  cumulative,  and  direct 
plots  to  be  discussed  in  the  next  section  may  be  output. 

If  one  introduces  a non- zero  noise  correlation 
matrix,  only  one  small  change  is  required  in  the  main 
program.  Due  to  core  storage  limitations  an  auxiliary 
program  is  used  to  store  on  tape  several  thousand  sets 
of  N numbers  (where  N is  the  number  of  stations  in  the 
simulated  network)  which  have  the  means,  standard 
deviations,  and  correlation  matrix  required  for  the 
noise.  Then  for  each  station  the  difference  is  taken 
between  the  randomly  generated  signal  value  and  the 
randomly  generated  noise  value,  read  from  tape.  If 
the  signal  is  greater  than  the  noise,  a detection  is 
declared . 


Simulation  Results 


First  let  us  discuss  a single  station  trying  to 
detect  all  P (or  LR)  signals.  Given  the  parameters 
in  Table  III,  an  answer  for  P may  be  converted  into 
an  answer  for  LR  simply  by  shifting  a "P"  graph 
5.80-3.47-0.33  units  (the  M -m^  threshold  difference 
at  LASA  for  Kuri 1- Kamchatka)  to  the  right.  Rows  3 and 
4 of  Table  V give  us  these  results,  and  they  are  de- 
rived from  Figure  8 as  indicated  in  the  third  column. 
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We  note  that  the  50 % incremental  value  is  in  exact 
agreement  with  the  mean  noise  levels  in  lable  III. 

This  is  as  predicted  in  the  theoretical  section  of  this 
report.  Also,  the  .08  m^  bias  in  the  asymptotic  por- 
tion of  the  observed  seismicity  curve  is  in  agreement 
with  theory.  In  fact  both  curves,  cumulative  and 
incremental,  are  in  complete  agreement  with  theory. 

Note  that  for  a single  station,  source  bias  cannot  be 
distinguished  from  other  sources  of  signal  variation 
and  that  the  correlation  of  noise  has  no  effect  on 
threshold.  We  note  that  the  9Ci  observed  incremental 
threshold  is  biased  by  0.22  with  respect  to  the  true 
magnitude  threshold  in  Table  IV  and  by  0.12  with 
respect  to  the  operational  magnitude  threshold. 


Next  we  consider  a single  station  which  must  first 
detect  on  P before  it  can  detect  on  LR  (we  note  that 
this  is  also  a good  model  for  two  separate  short- 
period  arravs,  e.g.,  NORSAR  and  LAS'X,  one  of  which  is 
prompted  by  the  other).  Figure  9 illustrates  the  LR 
thresholds  on  observed  M . For  the  parameters  given 
here,  comparing  rows  4 and  7 in  Table  V,  we  see  that 
there  is  very  little  effect  on  the  observed  Ms  thresh- 
olds when  the  P thre  .hold  lies  near  that  of  LR.  This 
mode  of  detection,  however,  offers  the  possibility  of 
other  kinds  of  plots,  such  as  those  shown  in  Figures  10 
and  11.  In  Figure  10  we  give  incremental  and  cumulative 
plots  of  the  number  of  events,  as  a function  of  observed 
m,  , for  which  LR  was  detected.  In  Figure  11 

D A 

the  fraction  of  events  for  given  observed  mb 


we  plot 
for  which 
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LR  was  detected.  Various  thresholds  in  terms  of 

observed  mb  for  detection  of  LR  may  be  obtained  from 

these  figures  and  are  given  in  row  6 of  Table  V, 

comparing  rows  6 and  7 we  see  that  the  m,  90%  thresh- 

b 

olds  arc  much  greater  than  the  thresholds,  although 
the  50%  points  arc  rather  close.  At  the  90%  level  the 
closest  agreement  of  all  may  be  found  between  the 
cumulative  mb  and  incremental  NETWORTH  thresholds  for 
true  magnitude  , 4.25  and  4.26  respectively.  This 
agreement  would,  of  course,  be  very  difficult  to  pre- 
dict from  physical  insight  (we  must  remember  that  if 

VV0-8  instcad  of  Ms=mb»  then  the  4.26  value  would 
become  3.46,  but  the  predictions  would  still  be  under- 
stood to  be  in  agreement). 

he  might  ask  what  would  be  the  effect  of  an  even 
Iom  P threshold  Comparison  of  rows  5 and  6 of 
Table  V shows  that  for  a P station  threshold  of  3.00 
instcad  of  3.47  there  is  almost  no  change  in  the  in- 
cremental 50%  and  90%  LR  thresholds  on  m.  . As  we 

b 

can  sec  by  comparing  Figures  12  and  10,  however,  the 
10 o incremental  threshold  would  be  significantly 
changed.  I t is  interesting  that  in  varying  the  P 
threshold  over  a large  range,  even  up  to  values  greater 
than  3.80  (the  LR  threshold),  the  single  station  direct 
probability  curve  in  Figure  11  did  not  change  notice- 
ably. This  result  is  reflected  in  Table  V.  It  would 
be  of  interest  to  prove  this  fact  theoretically.  We 
shall  sec  below  that  this  invariance  does  not  hold 
with  respect  to  network  thresholds. 
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Again,  these  s i n«,  1 e- s tat  i on  results  arc  unaffected 
hy  source  bias  or  hy  noise  correlation  between  stations. 
However,  there  woulJ  be  an  effect  if  there  were  corre- 
lation between  the  amplitude  levels  of  short- pc r i od 
and  long-period  noise;  this  topic  is  not  treated  in 
this  repo  i t . 

Next  we  move  to  the  subject  of  network  thresholds. 
The  thresholds  predicted  by  NITWORTH  have  already  been 
illustrated  in  figure  7.  We  investigate  first  the 
problem  of  four  stations  out  of  ten  detecting  without 
prompting.  I igure  13  is  the  appropriate  plot.  Although 
the  horizontal  axis  is  labeled  m.  , it  could  as  well  be 
M.  il  all  values  were  increased  by  0.33  (assuming 

W- 

We  note,  as  discussed  in  the  theoretical  section, 
that  there  is  a hump  in  the  seismicity-magnitude  curve 
in  contrast  to  the  single-station  case.  Hy  comparison 
with  figure  23,  for  which  no  source  bias  was  assumed, 
we  see  that  the  size  of  the  hump  is  reduced  by  the 
assumption  of  source  bias  (figures  22  and  23-31  corre- 
spond to  11  and  13-21  but  arc  without  source  bias; 

Table  VI  corresponds  to  Table  V).  We  illustrate  in 
figure  15  how  it  is  possible  to  estimate  erroneously 
the  seismicity  curve  and  to  obtain  a=1.2  instead  of 
0.8  in  the  seismicity-magnitude  relation.  Knowing 
that  the  hump  is  present,  we  can  draw  an  accurate 
seismicity  line  using  large  events  or  we  can  estimate 
the  slope  of  the  seismicity  lines  from  single-station 
incremental  plots  which  have  no  such  hump  and  thus 
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estimate  thresholds  which  are  reasonably  close  to  the 
predicted  ones.  We  find  that  the  seismicity  bias  is 
0.05  , approximately  what  equation  (8)  would  predict 

for  a single  station  with  a standard  deviation  of  the 
signal  equal  to  the  source  bias,  0.212.  This  is  a 
very  reasonable  result  in  that  for  large  magnitudes 
one  would  expect  the  average  observed  magnitudes  to 
be  very  close  to  the  true  magnitudes  except  for  source 
bias  (Figure  23  indicates  that  there  is  no  detectible 
seismicity-magnitude  bias  when  source  bias  is  absent). 

Comparison  of  rows  11  and  12  from  Table  V shows 
that  the  90%  incremental  threshold  is  in  good  agree- 
ment with  the  90%  true  magnitude  standard  NETWORTH 
threshold;  3.51  as  compared  to  3.52.  The  preferred 
comparison  would  be  with  NETWORTH  results  including 
source  bias,  for  either  true  or  operational  magnitide. 
i rom  Table  IV,  these  numbers  would  be  4.02  or  3.84 
minus  0.33  (3.69  or  3.51)  for  true  and  operational 
magnitude,  respectively,  instead  of  3.52;  thus  the 
agreement  is  quite  good  using  operational  magnitude, 
but  not  so  good  using  true  magnitude.  Note  that  the 
network  4/10  50%  point  is  in  almost  exact  agreement 
with  the  single  station  50%  point,  3.46  compared  to 
3.47.  This  must,  of  course,  be  regarded  as  a coinci- 
dence in  some  sense  since  the  50%  point  is  different 
for  1/10  or  2/10  station  detection. 

Finally  let  us  discuss  the  complete  network  system 
where  we  require  1/10  LR  detection  following  4/10  P 
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detection.  Figures  17a  and  27 , with  and  without  source 
bias  respectively,  give  the  incremental  and  cumulative 
M curves.  We  note  again  that  the  hump  is  noticeably 
reduced  by  the  presence  of  source  bias.  figures  16a 
and  26  show  that  plotting  the  incremental  and  cumula- 
tive curves  as  a function  of  the  m,  corresponding  to 
the  LR  detection  instead  as  of  a function  of  the  Mg 
value  greatly  reduces  the  hump  and  that  when  source 
bias  is  included  the  hump  apparently  disappear  com- 
pletely. The  direct  thresholds  are  determined  fiom 
Figure  11.  All  the  thresholds  calculatd  assuming 
source  bias  may  be  found  in  rows  15  and  16  of  Table  \. 
Rows  13  and  14  of  both  Table  V and  VI  were  calculted 
assuming  that  the  individual  station  P threshold  is 
lowered  from  3.47  to  3.00.  Comparison  shows  that  the 
90%  thresholds  are  hardly  changed  but  that  the  50% 
levels  change  by  as  much  as  0.21  magnitude  units. 

A 

Due  to  the  hump,  the  incremental  50%  and  90%  Mg 
thresholds  are  much  closer  together  in  the  simulation 
results  than  in  the  N'hTWORTII  results;  no  satisfactory 
fit  with  Table  IV  is  possible.  Also,  it  does  not  seem 
possible  to  fit  any  other  of  the  simulation  results 
closely  to  any  of  the  N FT WORTH  results.  Of  course  we 
should  emphasize  that  there  is  no  reason  to  expect  a 
particularly  good  fit;  the  whole  point  of  the  simula- 
tion is  to  derive  the  behavior  of  thresholds  obtained 
in  specified  empirical  ways  so  that  they  may  be  related 
to  the  true  or  operational  NUTWORT1I  thresholds.  In 
Table  V,  comparison  of  rows  15  and  16  to  rows  21  and  22 
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from  Figures  16b  and  17b  shows  that  if  one  has  p^-0.4 
in  addition  to  source  bias,  the  effect  is  almost 
unnot iceab le  at  the  SOI  level;  and  on  the  order  of 
0.05  magnitude  units  at  the  901  level. 

Similar  results  are  obtained  for  2/iC  and  4/10  LR 
detection,  and  these  results  are  given  in  Figures  18-21 
when  source  bias  is  included.  Compared  to  the  1/10 
thresholds,  the  2/10  ar.d  4/10  thresholds  increase  by 
0.10-0.15  with  each  step  increase  in  number  of  detecting 
stations  required.  This  is  in  rough  agreement  with, 
but  slightly  less  than,  the  NOT WORTH  differentials  of 
0.15-0.20. 

Figure  32  shows  the  probability  of  detection  of 
LR  by  1/10,  2/10  and  4/10  stations  as  a function  of 
the  observed  m^  determined  by  four  or  more  stations. 

The  curves  are  plotted  on  log-probability  paper  on 
which  the  cumulative  normal  probability  distribution 
function  would  plot  as  a straight  line.  The  failure 
of  the  2/10  and  4/10  curves  to  approximate  a straight 
line  suggests  that  fitting  observed  network  curves 
with  this  function  could  on  occasion  be  inappropriate. 
Note  that  this  function  fits  the  single-station  data 
well  though,  as  expected. 

The  main  point  which  the  reader  should  derive 
from  a detailed  study  of  the  foregoing  simulation  results 
for  a particular  ideal  network  is  that  any  casual 
approach  to  the  estimation  of  seismic  thresholds  for  a 
real  network  is  likely  to  be  in  error  by  several  tenths 
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real  n"aLnit"de  Unit'  ThC  b°St  “““  »f  action  when  a 
network  is  under  consideration  is  to  construct  a 

1° / y estlmatin«  its  parameters  from  the  most  accu- 

M hblG  d3ta  ^ t0  ' simulation  using 

MSBNET  on  this  model. 

Estimating  Relationships 

tion  oTm"!  ” !"“  th°  qUeSti°"  °f  thc  avaiua- 

V b relatl0ns'  Tabla  VII  gives  the  calculated 
onu  a lnterCepts  usi"S  regression  on  i,  regression 
“ maximum-I ihelihood  (Ericsson,  1971b)  for 
seven  different  types  of  m m . * ' 

1/10  7 / in  , ‘Vmb  data  cutoffs  and  for 

, ~/10,  4/10  LR  detection.  The  data  ,w 

aa  that  from  which  the  thresholds  pre!  t \ 

were  estimated,  and  the  simulation  values  for  th  0 e 

and  intercept  are  1 0 and  n n 10pe 

p i.u  and  o.o,  respectively 

r Hr 

rg„ttudi  estimates;  and  Column  HI,  using  data  )'£ 

• 5.0  as  low-magnitude  estimates.  Columns  IV  and  V 

encompass  the  magnitude  interval  from  6.0  down  to  that 
M above  whtch  lie  905  or  respectively,  of  the 

total  number  of  detected  events.  Column  VI  cuts  off 

" ded:::ti„Vih:;b  for  ^ r 

the  ClnTs sy™°tricai 

b S cut  offs  are  equal.  This  means  that  the 
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cutoff  lines  intorsert  A 

ntersect  at  the  mean  M -m  line.  For 

Columns  VT  I nn  rl  ....  s b 


l 


r i "c  ime.  hor 

Columns  VI,  and  VIM.  this  is  not  true.'  those 

columns  there  is  the  standard  all-inclusive  3.0-6  0 
an8°  tilC  however,  for  Column  VU  tl 

throVT0''  0n/'S  “ 8iV“  ^ th°  *.  ValU0  e<"*al  to 
the  Cl  " rCCt  CSh0ld  0n  V and  Column  VIII  by 

wh L*  “ "hiCh  lle  905  «f  events  for 


which  LR  was  detected 


nspection  of  this  table  shows  that  accurate  esti- 
sy  O"  °f  thc  Par“eters  a and  b are  obtained  using 

if  bothT  CUt°ffS  and  MXimun,-likelihood  estimation 
f both  low-magnitude  cutoffs  are  given  by  4.0  or  if  a 

u of  1S  made  at  the  ^ valuc  abov?  which  iio  a 

8Uo  Of  e\ ents  with  detected  i d rr  i 

erected  LR  (Columns  II.  iv  VI 

ccur  e resuUs  arc  als  „„  by  ^ J- 

tne  ro  iq  /- 1 . ~ r r * * c 


if  - ws  7 ‘CKicssion  on  M 


x the  Ms  value  above  which*  lie  001  of  alT  v t J h 
detected  LR  (Column  VII)  or  is  given  by  the  M value 
equal  to  the  LR  901  direct  threshold  (Column  VIII). 

accuraT  CaS°S  C°nSid°rcd  is  tegression  on  J. 

accurate.  However,  in  all  of  these  cases  there  is 

some  ias ; and  a simulation  approach  should  be  used  to 

Ptedtct  xt  and  thus  enable  one  to  remove  it. 

The  success  of  the  maximum-likelihood  approach  is 

and  r;1Sin8.bCCaUSC  aftCr  10"  c removed 

the  remaining  data  is  cut  off  so  as  to  retain  the 

ymme try  of  the  data  about  the  true  relation. 

then  the  maximum-likelihood  approach  77  be  expo  ted 
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to  he  unbiased  toward  either  too  large  or  too  small  a 
slope . 

Ihe  success  of  the  regression  on  Mg  is  perhaps 
surpiising.  It^certainly  has  nothing  to  do  with  greater 
reliability  on  Mg  values.  In  fact,  observed  Ms  values 
in  the  simulation  will  have  larger  variances  than  the 
observed  mb  values  when  LR  detection  is  not  as  good  as 
P detection.  However,  such  a cutoff  will  ensure  that 
at  each  value,  in^  values  will  be  scatter  approxi- 
mately normally  about  the  mean  line  (this  is  not  exactly 
true,  since  we  have  seen  in  the  theoretical  section 
that  realistic  seismicity  shifts  the  mean  of  the  normal 
distribution).  On  the  other  hand,  the  same  cutoff  will 
ensure  that  near  this  cutoff  the  Mg  values  are  quite 
dramatically  not  normally  distributed  about  the  mean 
line,  for  fixed  mb . And,  of  course,  such  a normal 
distribution  is  a requirement  for  validity  of  regres- 
sion calculations. 

The  point  of  this  section  is  that  one  should  simu- 
late the  nctuork  under  consideration  and  so  discover  by 
trial  and  error  which  techniques  for  estimation  of  the 
Ms " mh  line  result  in  the  smallest  bias.  Then  these 
techniques  should  be  applied  to  the  real  data,  and 
any  expected  errors  can  he  allowed  for  in  the  inter- 
pretation. Results  from  our  particular  network  siru- 
lation  indicate  that  maximum- 1 ikcl ihood  estimation  with 
symmetrical  cutoffs  on  the  data  at  the  90$.  LR  threshold 
is  an  accurate  technique  and  that  regression  of  m on 

Ms  Wlth  a 9()°6  threshold  cutoff  on  Mg  only  is  a good 
alternative  technique. 
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In  an  earlier  theoretical  section,  a method  of 
estimating  the  true  M.-m^  line  was  outlined;  essentially 
it  involves  the  use  of  equation  (25).  This  method  is 
.alid  for  single-station  plots  only,  and  it  also  re- 
quires more  precise  knowledge  of  detection  probabilities 
than  maximum  likelihood  and  standard  regression  methods 
do.  An  illustration  of  the  method  is  worthwhile  though. 
Using  MSBNIiT  simulated  data  for  one  station  of  t’  e pre- 
viously employed  network,  an  incremental  plot  of  number 
of  events  versus  m^ , shown  in  Figure  33,  was  made  to 
determine  the  seismicity-magnitude  relation  and  thus 

the  expected  number  of  events  N.  at  each  m,  . (The 

i b . v 

low  P threshold  case  with  ^=3.00  was  used.}  The 
surface-wave  and  body-wave  detection  probabilities  for 
this  station,  relative  to  station  magnitude,  are  given 
by  equation  (2)  with  on=0.2  and  50%  probabilities  at 
3.80  and  3.00,  respectively.  In  practice,  these  prob- 
abilities could  be  estimated  from  noise  values  or 
incremental  detection  plots.  This  information  was 
applied  in  the  evaluation  of  equation  (25).  A value 
of  e=0.1  was  selected  arbitrarily  for  one  set  of 
evaluations . Also  a value  of  e=0.35,  as  determined 
by  equation  (23),  was  used  to  attempt  to  estimate 
exactly  the  true  M -in.  line.  Figure  34  shows  the 
MSBNLT  data  and  the  calculated  (Mg  ,m^  ) points.  No 
evaluations  were  made  for  mb<3.3,  whicft  is  roughly 
the  .95  detection  probability  point  for  body  waves  as 
seen  in  Figure  33;  and  no  summations  were  made  which 
included  a surface-wave  detection  probability  less 
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than  .16,  at  Ms=3.6.  The  results  appear  to  be  excel- 
lent, with  the  e=0.35  points  closely  defining  the  true 
Mg-HU  line.  The  line  has  been  estimated  down  to 
= 3.3,  the  95"o  body-wave  threshold  for  the  station. 
Note  that  the  method  should  work  reasonably  well  even 
if  the  true  Ms*m^  relation  were  not  linear.  In  that 
case  the  seismicity-magnitude  relation  would  not  be 

A A 

log- linear  and  the  density  function  f(Ms|m,)  given  by 
equation  (20J  would  no  longer  be  strictly  correct; 
but  the  deviations  should  not  be  large  enough  to 
nullify  application  of  the  procedure,  illustrated  here 
for  a linear  M^-m^  relation,  to  somewhat  non-linear 
relations . 

A closely  related  M -m^  estimation  procedure 
which  does  not  require  estimation  of  the  seismicity- 
magnitude  relation  nor  station  M thresholds  is  pos- 
sible for  mb  values  above  the  90?o  detection  threshold. 

In  this  case  one  s imply  sums  down  to  an  M which 

s . 

gives  a fixed  fraction  c of  the  events  actually  ob- 
served at  m^  . This  assumes  perfect  detection  down 
to  Mg  . The 1 resulting  curve  fit  to  the(Ms  ,mb  ) pairs 
may  tl\en  be  shifted  down  by  the  di  f fcrence1at  ftigh 

A 

magnitudes  between  the  M.  's  obtained  using  e and  the 

a ^ i 

M 's  obtained  using  an  e corresponding  to  the  true 
s"mb  (cm  ~ 0.5  might  not  be  too  inaccurate  for 

most  practical  purposes! . We  have  evaluated  this  pro- 
cedure using  program  MSBNliT  and  find  that  it  is  appa- 
rently accurate,  although  somewhat  unstable,  for  the 
case  c*0.1,  which  was  chosen  so  that  the  estimation  of 
the  line  could  be  extended  to  low  magnitudes. 
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Comparison  of  Simulation  with  LASA  M Data  for  the 

Kamchatka- Kuri 1 Region 

Figure  35  shows  the  incremental  and  cumulative 
curves  for  the  SAAC-LASA  daily  summary  for  geographic 
regions  217-222  in  the  Kamchatka- Kuri 1 area  for 
approximately  six  months  in  1972.  The  a value  is 
approximately  0.8  as  in  the  simulation.  The  incremen- 
tal and  cumulative  thresholds  from  this  figure  are 
given  in  Table  VIII.  Mack  and  Robertson  (1973)  gave 
a cumulative  curve  from  which  it  is  possible  to  obtain 
the  cumulative  501  and  901  points  for  LASA  LR  detection 
in  the  Kamchatka-Kuril  region  for  events  having  a P 
detection  as  reported  in  the  LASA  bulletins.  They 
also  showed  a curve  of  direct  probability  of  LR  detec- 
tion as  a function  of  LASA  m^. 

In  choosing  our  simulation  parameters,  we  had 
these  LASA  data  in  mind;  therefore  we  adjusted  the 
mean  s ingle- station  short-period  and  long-period  noise 
levels  so  that  the  difference  in  simulated  thresholds 
would  be  as  observed  at  LASA,  and  we  also  chose  the 
absolute  levels  of  the  noise  values  so  that  the  simu- 
lated thresholds  were  in  agreement  with  the  observed 
/alues  (to  within  .01  magnitude  unit).  Finally, 
setting  M am^-0 . 8 instead  of  Ms=mb  as  in  the  simulation 
enables  us  to  present  in  Table  VIII  a comparison 
between  observed  and  simulated  threshohls.  In  general 
we  feel  that  the  agreement  is  within  the  experimental 
error.  More  conclusive  tests  of  the  agreement  of 
simulations  with  observation  shou’d  be  possible  using 
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the  LI\.  data.  We  shall  treat  this  problem  in  future 
reports.  It  \%ill  he  important  to  see  if  network  in- 
cremental detection  curves  for  small  source  regions 
do  indeed  have  humps. 


REFINEMENTS  IN  CALCULATED  THRESHOLDS 


S/N  Ratio  Required  for  Detection 
A variation  of  equation  (1) 


I’dU) 


logA  - 


2 

(°“  + 


u + logr 


(28) 


is  most  often  used  in  predicting  thresholds;  here  r is 
the  S/N  ratio  required  for  detection.  An  assigned  value 
of  r must  depend  on  the  frequency  band  of  the  signal 
and  of  the  noise,  and  on  the  analysis  procedure.  If 
events  having  nearly  the  same  frequency  content  as  the 
background  noise  are  being  searched  for  without  know- 
ledge of  their  time  and  frequency  of  occurrence,  then 
a rather  high  r may  he  set  for  calculating  detection 
thresho I . On  the  other  hand,  if  events  with  frequency 
content  different  from  the  recorded  noise  and  with 
accurate  predicted  arrival  times  are  being  search  for, 
a rather  low  r must  apply.  It  is  mandatory  that  r be 
determined  via  a controlled  experiment  when  calculated 
detection  probabilities  arc  to  be  compared  with  empiri- 
cal probabilities  determined  in  analysis  of  new  seismic 
data.  Such  an  experiment  is  reported  by  Geotech  (1966) 
f o ’ short-period  detection  at  the  VELA  observatories 
uicrc  a S/N  ratio  (peak- to-peak  signal  over  peak-to- 
peak  noise)  of  roughly  1.5  seemed  appropriate  for  detec- 
tion at  the  501,  level.  Another  experiment  is  reported 
by  von  Seggern  (1974)  for  the  Long-Period  Experimental 
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network  where  a S/N  ratio  (peak- to-peak  signal  over 
peak-to-peak  noise)  of  roughly  1.0  seem  appropriate 
for  LR  detection  at  the  50 % level.  Note  that  the 
choice  of  definition  of  y,  pcak-to-pcak  or  rms,  affects 
the  value  of  r.  There  is  approximately  a factor  of 
six  difference  in  the  two  definitions,  since  if  the 
noise  is  Gaussian,  a peak-to-pcak  noise  excursion  six 
times  greater  than  the  rms  is  expected  once  for  roughly 
every  7 minutes  of  long-period  data  sampled  at  1 pt/scc, 
or  seven  times  greater  thin  the  rms  once  for  roughly 
every  33  minutes. 

Background  Noise  Level 

The  value  of  y in  equation  (28)  should  be  deter- 
mined empirically  from  actual  noise  recorded  at  the 
site.  It  is  often  taken  as  the  rms  value  of  a noise 
sample  or  average  rms  of  many  noise  samples  (c.g., 

Benno,  1972).  For  the  purpose  of  detecting  signals 
however,  there  is  only  a tenuous  connection  between 
the  rms  value  and  the  noise  peaks  which  effectively 
set  the  detection  threshold.  The  rule  of  thumb  for 
multiplying  rm ; by  six  to  get  peak-to-peak  noise  maxima 
only  applies  if  the  noise  is  Gaussian,  this  requirement 
being  especially  important  for  extreme  values  of  the 
distribution.  Therefore  a simple  visual  assessment  of 
peak  noise  values  (Gcotcch,  1966j  von  Scggcrn,  1974) 
is  at  least  as  good  as  and  in  some  cases  undoubtedly 
better  than  rms  measurements.  This  approach  estimates 
exactly  what  maximum  noise  amplitudes  the  human  or 
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raac'tine  detector  is  up  against  to  reliably  pick  signals 
without  creating  an  unacceptable  false  alar.ii  rate. 


Another  matter  to  be  considered  in  determining 
p in  equation  (28)  is  the  presence  of  interfering 
events  on  seismic  recordings.  Many  events  which  prob- 
ably would  have  been  detected  in  normal  background 
are  masked  by  larger  signals  from  another  event. 

These  interfering  events  should  logically  be  considered 
as  the  upper  values  of  noise.  If  p is  determined 
solely  from  quiescent  periods,  the  resulting  predic- 
tions of  station  or  network  capability  using  (28)  will 
necessarily  be  lower  than  actually  achieved.  This 
bias  is  not  insignificant.  In  Figure  36  we  show  two 
noise  distributions  for  the  Long-Period  Experimental 
station  at  Charters  Towers,  Australia.  Both  distribu- 
tions represent  the  picking  of  maximum  peak-to-peak 
excursions  of  "noise"  with  period  near  20  seconds  on  the 
vertical  component  in  a 40-minute  window  at  roughly 
equally-spaced  dates  in  1972;  the  samples  with  the 
lower  mean  were  taken  always  at  seismically  inactive 
times,  and  the  samples  with  the  higher  mean  were  taken 
at  a fixed  time  of  day  no  matter  what  was  present  on 
the  recording.  The  logarithm  of  the  means  differ  by 
nearly  .3  magnitude  unit.  We  suggest  that  the  higher 
p is  the  appropriate  one  to  use  in  predicting  thresholds 
if  interfering  events  cannot  be  eliminated  in  detection 
studies . 


SUMMARY  AND  CONCLUSIONS 


In  order  to  properly  define  seismic  thresholds, 
it  was  necessary  to  define  three  levels  of  seismic 
magnitude:  true,  ope  rational , and  observed.  The 

observed  magnitudes  arc  continually  at  our  disposal, 
but  we  must  determine  thresholds  on  operational  magni- 
tude for  the  purpose  of  fair  comparisons  of  detection 
capability  among  stations  and  networks.  A further 
refinement  into  the  frame  of  true  magnitude  ties 
thresholds  to  some  absolute  indicator  of  event  magni- 
tude, such  as  total  energy  release  in  explosions  or 
earthquakes.  Most  of  the  difficulties  in  treating 
thresholds  can  be  traced  to  the  fact  that  ratural 
seismicity  results  in  more  events  scattering  into  a 
certain  magnitude  from  lower  magnitudes  than  from 
higher  ones  end  to  the  fact  that  there  is  increasing 
magnitude  bias  with  decreasing  magnitude. 

Many  of  our  problems  in  threshold  determination 
can  be  handled  analytically.  We  established  in  this 
manner  ^nat: 

1.  for  a single  station,  the  50%  incremental 
threshold  magnitude  using  observed  magnitudes  is  the 
same  as  the  50%  direct  threshold  determination  using 
operational  magnitudes;  but  the  90%  and  10%  arc  lower 
and  higher,  respectively. 

2.  The  single-station  incremental  curve  on  observed 
m smoothly  dies  off,  but  the  network  curve  has  a hump 
near  the  threshold. 
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3.  I he  sc i smici ty-magni tude  relation  determined 
using  observed  magnitude  will  always  be  displaced 
upward  in  magnitude  compared  to  the  true  seismicity. 

4.  Cumulative  threshold  magnitudes,  either  net- 
work 01  s ingle-s  t r.t  i on  t will  always  be  lower  than  the 
incremental  ones. 

5.  1 he  distribution  of  magnitudes  aiong  a straight 
line  (horizontal  or  vertical)  through  a population  of 
^s"mb  ^ata  i-s  dependent  on  the  seismicity-magnitude 

l elation  and,  although  normal,  aas  varying  parameters. 

»\e  have  presented  the  method  for  converting  seismic 
thresholds  for  LR  from  to  scales  and  vice  versa. 
Because  it  is  necessary  to  invoke  several  estimated 
parameters  in  this  procedure,  such  as  the  amount  of 
source  bias,  the  extent  of  path  effects  on  amplitude, 
the  slope  and  intercept  of  the  M^-m^  relation,  and 
r-he  slope  and  intercept  of  the  seismicity-magnitude 
relation,  we  feel  this  conversion  should  be  done  with 
care,  lest  results  stray  from  predictions. 

he  have  assessed  the  effect  of  source  bias  and 
noise  correlation  on  detection  thresholds  using  a 
modified  version  (M  IWCORR)  of  our  standard  program 
for  predicting  network  detection  capability.  Results 
show  that  network  capabilities  arc  altered  insignifi- 
cantly by  noise  correlation  and  that  predicted 
thresholds  for  true  and  operational  magnitudes  typically 
differ  by  roughly  .1  magnitude  unit. 

An  accurate  method  for  reconstructing  the  true  M,-m^ 
line  1 rom  single  station  data  was  derived  theoretically 
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using  the  statistical  properties  of  the  observed  data 
compared  to  the  total  population  of  possible  data. 

A simulation  computer  program  was  written  to  gain 
insights  in  such  matters  as  network  threshold  problems 
not  analytically  tractable,  effects  of  correlated  sig- 
nals and  noise,  feedback  between  P and  LR  thresholds, 
and  estimation  of  M -in^  relationships.  Through  the 
analysis  of  an  ideal  network  we  have  attempted  to 
illustrate  the  range  of  situations  which  arise  in 
threshold  determination.  Comparison  of  program  MSB NET 
results  with  those  of  real  networks  enables  one  to 
infer  certain  parameters  for  the  network  and  to  make 
threshold  corrections  to  the  empirical  data. 

Simulation  results  agreed  with  theoretical  pre- 
dictions wherever  compar;son  was  possible.  The  simu- 
lation experiment  established  these  important  facts: 

1.  Varying  the  P threshold  iocs  not  affect  the 
single  station  direct  LR  threshold  on  m^  but  does 
affect  the  network  threshold.  Incremental  and  cumula- 
tive thresholds  can  be  strongly  affected  in  either 
case . 

2.  The  hump  in  the  network  LR  incremental  detec- 
tion curve  is  reduced  by  assuming  significant  source 

.as  and  also  by  plotting  versus  m^  rather  than  . 

3.  The  90 % s ingle- stat ion  LR  thresholds  in  terms 
of  in^  are  typically  much  greater  than  the  thresholds 

in  terms  of  M t assuming  M.  = in,  ) oi  the  order  of  0.4 

S'  s a 

magnitude  units.  In  general,  really  substantial  effects 
are  possible  and  care  must  be  exercised  when  comparing 
thresholds  determined  in  different  ways. 
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4.  Noise  correlated  between  stations  at  realistic 
levels  (p-0.4)  has  an  almost  negligible  effect  on 
thresholds  in  the  presence  of  realistic  levels  of 
source  bias. 

5.  1 he  seismicity  bias  for  large  magnitudes  is 
negligible  if  network  magnitudes  are  used,  except  in 
the  presence  of  source  bias  in  which  case  it  may  be 
on  the  order  of  0.1  magnitude  unit. 

6.  lor  realistic  network  parameters,  compared  to 
the  1/10  thresholds,  the  2/10  thresholds  increase  by 
0.10-0.15  magnitude  unit  and,  compared  to  the  2/10 
thresholds,  the  4/10  thresholds  increase  by  a similar 
amount.  This  is  in  rough  agreement  with,  but  slightly 
less  than,  the  NET WORTH  differentials  of  0.15-0.20. 

7.  Use  of  the  error  function  as  a probability  of 
detection  curve  for  a network  can  in  some  cases  be  a 
poor  approximation  to  the  curve  deduced  by  simulation. 

8.  The  true  Ms - m^  relationship  is  accurately 
estimated  by  the  maximum  likelihood  method  if  a suit- 
able low  cutoff  is  made  on  the  and  m^  values. 
Regression  of  m^  on  also  yields  accurate  results 

if  that  10%  of  tne  events  with  the  lowest  M values 

s 

are  excluded  from  the  calculation. 

As  an  example  of  the  usefulness  of  a simulation 
approach  to  threshold  problems,  good  agreement  was 
achieved  between  simulated  results  from  MSBNET  and 
LASA  data  for  LR  and  P detection. 
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In  regard  to  predicting  seismic  thresholds,  we 
reiterate  the  importance  of  determining  realistic 
background  "noise"  values  which  take  into  accourt 
seismic  waves  from  known  sources  and  also  the  impor- 
tance of  determining  the  S/N  ratio  required  for  detec- 
tion in  a specified  context. 
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Simulation  Parameters  for  Network 


TABU:  IV 

Predicted  Thresholds  for  the  Simulated  Network 
Looking  at  Kuri 1-Kamchnt  a Region 
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Note:  See  Table  III  for  network  parameters. 
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Predicted  and  Observed  LASA  Kuril-Kamchatka  Thresholds 
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I igurc  1.  Incremental  f rcqucncy- magni tudc  plots  for  a 
hypothetical  data  set  with  true,  s tation- obsc rvcd , and 
network-observed  magnitudes  as  the  basis  (no  source  bias) 
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Figure  2,  Typical  plot  for  earthquakes. 
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Figure  5.  Theoretical  effect  of  correlated  signals 
noise  on  a two-station  network. 


Figure  t>.  Predicted  effect  of  correlated  noise  on  the 
Long  Period  experimental  network. 


Figure  7.  NETWORTH  capability  estimates  for  10  LASA-type 
stations  equally  spaced  about  Kurii-Kamchatka  with  no 
source  bias  present. 


NUMBER  OF  EVENTS 


Figure  11.  Direct  LR  thresholds  for  the  simulation 
network  given  that  >4  station  P detection  has  been 
declared  - source  bTas  present. 
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Figure  13.  Network  P (or  LR)  thresholds  on  m,  (or  M ) 
for  >4  out  of  10  detecting  with  source  bias. 
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f igure  14.  Network  LR  thresholds  on  m,  for  >_1  out  of 
10  detecting  with  source  bias  and  low  TP  threshold  (3.00) 
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figure  16b.  Network  LR  thresholds  on  m,  for  >1  out  of  10 
detection  with  source  bias  and  standard0?  threshold  (3.47) 
and  correlatbd  noise  (p  =.4), 
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figure  17b.  Network  LR  thresholds  on  M for  >1  out  of  10 
detecting  with  source  bias  and  standardsP  threshold  f3  471 
and  correlated  noise  (pn=.4).  ' 
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f igure  18.  Network  LR  thresholds  on  m,  for  >2  out  of  10 
detecting  with  source  bias  and  standard  P thTeshold  (3.47). 
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R thresholds  on  m,  for  >4  out  of  10 
bias  and  standara  P threshold  (3.47). 
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f igure  21  „ Network  LR  thresholds  on  M for  >4  out  of  10 
detecting  with  source  bias  and  standard  P threshold  (3.47). 


-112- 


Figure  22.  Direct  LR  thresholds  for  the  simulation  network 
given  that  a >4  station  P detection  has  been  declared  - no 
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Figure  23.  Network  P (or  LR)  thresholds  on  m,  (or  M ) for 
>4  out  of  10  detecting  without  source  bias. 
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rWnrt-25*  •N,etwork  LR  thresholds  on  M for  >1  out  of  10 
detecting  without  source  bias  and  withslow  P'threshold^.OO) 
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Figure  26.  Network  LR  thresholds 
detecting  without  source  bias  and 
threshold  (3.47). 
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Figure  30.  Network  LR  thresholds 
10  detecting  without  source  bias 
threshold  (3.47). 
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rigure  31.  Network  LR  thresholds  on  M for  >4  out 
threshold' (3.4 / K°Ut  S°UrCe  bUs  and  with  5tandard 
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Figure  32.  Network  LR  probability  of  detection  curves 
o n ni|^  & 
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LOGARITHM  OF  NOISE  AMPLITUDE  IN  mp  ZERO  TO-PEAK 

Figure  36.  Distribution  of  CTA  long-period  noise  amplitude 
picked  on  the  high-gain  vertical  component  in  1972. 


APPENDIX  I 

DERIVATION  OF  THE  PROBABILITY  DENSITY  OF  OBSERVED 
MAGNITUDE  FOR  A TWO- STATION  NETWORK 
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appendix  I 


Derivation  of  the  Probability  Density  of  Observed 
Magnitude  for  a Two-Station  Network 

When  two  stations  detect  events  of  fixed  opera- 
tional magnitudes  m,  the  probability  density  of  the 
combined  observed  magnitudes  is  given  by  the  convolu- 
tion of  the  individual  densities  if,  as  is  clear  from 
physical  considerations,  the  individual  densities  are 

uncorrelated: 


f*  * (n)  ' f f;  (m  -mjf*  (m  )dm 
m1+m2v  cJ  * m2  c I m^  i i 

where  m^m^m^  We  will  convert  to  average  magnitude 


(Al) 


in  later. 

The  expression  for  f(m)  for  fixed  m at  a single 
station  has  been  given  by  Herrin  and  Tucker  (1972)  as 


f Cm) 


1 


exp 


. (A2) 


Substituting  (A2)  in  (Al)  and  rearranging  the  exponents 
results  in  : 


Preceding  page  blank 
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we  express  the  integral  in  (A3)  as: 
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Y/t  - z 
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( 2 TT ) 


T77 


~ z 2 / 2 . 
e dz , 


(A4) 


Using  the  relation  l-4>(-x)  = 4>(x),  we  can  separate 
the  integral  into  two  parts: 


z + y/t1 

1 c"z2/2d-  f t 
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(271) 
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(2tt) 


T77 


(AS) 


A theorem  from  Iillison  (1964)  provides  the  evaluation 
of  the  first  integral  in  (A5) : 


z + y/t 
' 1/x 


r Y 

® — 2 177 

(t  + l)1/z 


(A6) 


Another  theorem  from  Zachs  and  liven  (1966)  provides 
the  evaluation  of  an  integral  similar  to  the  second 
one  in  (A5) , except  that  their  means  in  the  cumulative 
normal  distributions  are  of  the  same  sign  as  well  as 
being  equal;  an  obvious  substitution  though  of  a 
positive  for  a negative  sign  can  be  carried  through 
their  derivation  to  produce: 
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( A7) 


where  4>(-p,  p , p)  is  the  bivariate  normal  distribution 
with  means  -p  and  p and  correlation  coefficient  p. 

Expanding  the  bivariate  normal  distribution  in 
terms  of  derivatives  of  the  normal  density  function 
gives: 
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The  product  of  even-ordered  derivatives  in  the  above 
series  is  positive  while  that  of  odd-ordered  derivatives 
is  negative,  so  that: 
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Now  using  equations  A9 , A7,  and  A6  in  A4,  we  can  ex* 


This  expression  goes  asymptotically  to  0 as  y-*-°° 
and  asymptotically  to  1 as  y-*°°;  its  value  at  y=0  and 
t*l  is  The  series  summation  is  rapidly  convergent 
for  all  y and  all  t.  Expressing  y and  t in  terms  of 
the  original  variables  and  then  substituting  (A10)  in 
(A3)  gives  the  density  function  for  mc: 
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(All) 


'«'^n 
probability  axiom: 


where  a "= (o^+o^/ 2) . I f we  define  m=mc/2,  then  by 


f (m)  * 2f*  * (2m). 

v m^  + n^ 


Applying  this  to  (All)  gives  the  desired  probability 
density  function  for  observed  average  magnitude  when 
both  stations  of  a two-station  network  detect: 
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APPENDIX  II 


New  NETWORTH  Capabilities  for  Handling  Source 
Bias  and  Noise  Correlation 

A.  IDENTIFICATION 

Title:  NETWCORR 

Programmer:  Modification  by  H.  Husted  of  NETWORTH  by 

Mark  Wirth 

Date:  20  September  1973 

B.  PURPOSE 

To  compute  the  network  magnitude  bias  for  the 
threshold  magnitudes  computed  by  the  standard  NETWORTH 
program.  Also  to  evaluate  the  effect  on  thresholds  of 
noise  and  signal  amplitudes  which  are  correlated  be- 
tween stations. 

C.  USAGE 

1.  This  is  a FORTRAN  IV  Program  for  the  360/44.  The 
principle  program  writeup  is  that  of  Wirth  (1970)  and 
having  it  at  hand  is  necessary  in  order  to  use  this 
one . 

2.  Modifications  to  NETWORTH  control  card  formats: 

If  bias  calculations  are  desired,  a 1 is  placed 
in  column  80  of  the  A card.  If  in  addition,  the  run 
is  to  be  made  with  correlated  signal  and/or  noise;  a 
2 is  placed  in  column  80.  Bias  will  be  computed  in 
this  case  also. 


Preceding  pege  blank 
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If  there  is  a 1 or  2 in  column  80  of  the  A card, 
then  immediately  following  must  come  a card  with  (NSD, 
NRM,  NRANDS , NREP,  NS,  NN,  NREPK)  (7110).  NSD  is  the 
number  of  stations  which  must  detect  before  a magnitude 
is  calculated.  If  correlated  signal  and/or  noise  are 
to  be  considered,  this  is  reset  internally  to  the  num- 
ber of  stations  required  for  detection,  found  in  columns 
9 and  10  of  the  A card.  NRM  is  the  number  of  randomly 
generated  event  threshold  magnitudes  which  must  be 
calculated  before  the  average  bias  is  computed.  NRANDS 
is  a random  number  to  start  the  random  number  generator. 
Default  lo  77773777.  These  are  the  only  additional 
parameters  needed  for  the  case  of  a simple  bias  calcu- 
lation. 

Tf  one  wishes  also  to  evaluate  the  effect  of  cor- 
related noise  and/or  signal;  then  the  meaning  of  NRM 
changes  slightly  and  the  additional  parameters  NREP, 

NS,  NN  and  NREPK,  together  with  covariance  or  correla- 
tion matrices  of  the  signal  and  noise  must  be  road  ir». 

NREP  is  . .amber  of  events  which  are  randomly 
generated  at  each  epicenter  in  an  attempt  to  have  NRM 
of  them  detected  so  that  a reliable  probability  of 
detection  e t.imate  may  be  obtained.  NkEPK  is  the 
maximum  number  of  times  to  repeat  NREP  trials.  (A 
limit  must  be  set  since,  for  low  magnitudes,  detections 
may  not  occur  even  after  large  amounts  of  computer  time 
have  been  used.  Some  care  in  choosing  NREP  and  NREPK 
is  needed  to  minimize  the  probability  biases  which 
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would  result  from  cutting  off  the  calculations  just 
when  NRM  detections  have  been  achieved.)  NS  and  NN 
are  0,  1,  or  2 as  the  signal  and  noise  are:  a)  uncor- 

related; b)  correlated  and  covariance  matrix  is  to  be 
read  in  ; or  c)  correlated  and  correlation  matrix  is 
to  be  read  in.  The  lower  triangular  portion  of  the 
covariance  or  correlation  matrices  are  read  in,  signal 
matrix  first.  (See  the  attached  data  forms.)  If  NS 
or  NN=0,  the  corresponding  diagonal  matrix  must  not  be 
read  in. 

3.  Timing 

The  program,  especially  with  correlated  noise  and 
signal,  can  be  very  time-consuming.  Restriction  of 
magnitude  range  and  careful  selection  of  epicenters, 
together  with  preliminary  test  runs  is  recommended. 

4.  In  general  all  options  of  NETWORTH  and  NETPLOT 
are  retained  with  the  exception  that  one  cannot  do 
subsets  with  correlated  signal  and/or  noise.  A default 
to  the  standard  NETWORTH  pregram  occurs  if  all  new 
parameters  are  zero. 

D.  METHOD 

If  a single  bias  calculation  is  desired  a magni- 
tude differential  AmjJ;-*  is  generated  for  station  i, 
epicenter  j. 

Am^  = N(0  ,a?s) 

where  a.  is  the  standard  deviation  of  the  signal  at 
is 
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station  i.  Then  if  is  the  amplitude  for  the 

threshold  magnitude  for  station  i,  epicenter  ; and  vm 
is  the  noise  level  we  compute 


*(xJJ)  where 


where  $ is  the  cumulative  normal  probability  distribu- 
tion, o.  is  the  noise  standard  deviation,  and  r is  the 
in 

signal- to-noise  ratio  required  for  detection. 

Then  a second  random  number  U,  , uniform  on  the 

^ i i 

interval  0<U^,<1  is  generated,  and  if  £ we  detect 

and  add  Am*^  to  the  sum  over  detecting  station: 

!Vj. 

Afte  ' all  stations  have  been  considered,  we  divide  by 
the  number  of  stations  detecting  if  ^NSD  yielding  A^ . 
Repeat  until  NSD  or  more  stations  have  detected 
NRM  times  and  save  the  average  (l/NRM)^Am^  as  the  bias 
for  epicenter  j„ 

If  noise  and/or  signal  correlation  exists  then 
for  each  "random  event"  sets  of  random  numbers  with  the 
required  mean  and  covariance  structure  are  generated  by 
the  methods  outlined  by  Shumway  and  Blandford  (1970). 

Then  detection  occurs  at  each  station  if  {logA+Am- (y+logr) } 
is  greater  than  zero.  Sufficient  iterations  may  be 
taken,  as  outlined  in  C-2  above  to  ensure  stable  esti- 
mates of  probability  of  detection  and  bias. 
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F.  DATA  CARD  FORMATS. 

Attached. 


DATA  FORM  netwcorr  [""page  i I name 


ESO  A77  2/64 


HHODMinN  WHOd  viva 


ESO  A77  2/64 
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APPENDIX  III 

Writeup  for  the  Program  MSBNET 

1.  Program  Name  - MSBNLT 

2.  Purpose  - Analysis  of  detection  statistics  of 
seismic  networks  - see  reference. 

3.  Procedure: 

livery  network  simulation  begins  with  an  initial 
set  of  events,  each  specified  by  its  "true"  magnitude 
Nb.  The  range  of  body  wave  magnitudes  under  considera- 
tion increases  from  a minimum  of  Ml  (Input)  in  incre- 
ments of  0.1  magnitude  unit,  and  the  expected  number 
of  events  in  an  interval  [m^m^]  of  these  "true"  values 
is  given  by: 

EIN(M))  • f2  Ao“("14»d.  (Ij 

ml 

A , u > 0 

where  A and  a are  input  parameters.  The  total  number 
of  events,  A/~  must  not  exceed  5000.  The  Nb  are 
generated  by  seeing  if  a random  number,  uniform  on  the 
interval  (0,1)  is  greater  than  the  total  number  of 
expected  events  up  to  Nb  divided  by  M-p,  the  total 
number  expected,  and  less  than  the  expected  number  up 
to  M.  , divided  by  M~.  If  so,  one  event  of  magnitude 
Nb  is  generated.  Once  the  total  number  of  expected 
events  has  been  generated,  the  program  proceeds.  Of 


-All  1-2- 


course  some  magnitude  cells  will  be  overfilled  and  some 
underfilled,  just  as  with  real  seismicity  data. 

Once  the  sample  seismicity  has  been  determined, 
the  following  steps  are  carried  out  for  each  station: 

1.  Scatter  (due  to  variations  in  source  regions, 
transmission  paths,  etc.)  is  simulated  by  generating 
a randomized  population  of  station  amplitudes  from  the 
initial  magnitude  set  as  follows: 


am. 

1 


+ b ♦ y,  (0, 


a2  ) 
ss; 


+ e. 


CO,  ) 
. * ms' 


A 


m 


b 


k 


mi  * 6k(0-  °sb)  * nkt0>  0mb> 

i = 1,  2,  ...,  100 

k * 1,  2,  ...,  N(M.) 


(2) 


where  a is  the  slope  and  b is  the  intercept  of  the 
true  Ms-mb  relation  and  the  yk,  6k,  and  nk  are 
random  normal  variates  with  the  indicated  means  and 
variances.  The  standard  deviations  o b and  o may 
vary  as  input  parameters  from  station  to  station,  while 
°ms  and  °mb  are  constant,  y and  6 are  different  for 
every  event,  but  the  same  for  every  station.  That  is, 
e and  n are  the  "source  bias"  and  y and  6 are  the 
"path-station  bias". 
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2.  Least-squares  regression  lines  on  mb  and  ftg 
fitted  to  the  values  obtained  from  equation  (2).  In 
addition,  a maximum  likelihood  linear  fit  (see  writeup 

A » 

of  program  MAXLIK)  is  computed  for  M * .n^  which  allows 
for  the  presence  of  normally  distributed  errors  in 
both  variables.  The  slopes  and  intercepts  are  printed 
and  also  accumulated  to  provide  averages  over  all 
stations. 

3.  The  population  is  searched  to  find  the  number  of 

A 

events  scattered  into  each  M magnitude  interval: 

N : (aMi  - .05)  < ftR  < (aAL  + .05) 

i = 1,  2,  ...,  100 

A 

then  again  on  nu  : 

Nb  : (M.  - .05)  < mb  < (M.  + .05) 

1 1 i = 1,  2,  ...,  100 

These  numbers  are  printed,  stored  individually,  and 
also  accumulated  over  all  stations. 

4.  To  provide  a detailed  measure  of  the  shape  of  the 

A 

scattered  population,  the  average  mb  is  computed  for 
each  fixed  value  of  ft  . These  averages  a-e  printed 
and  accumulated.  In  addition,  the  average  M for 
fixed  mb  is  computed,  printed,  and  accumulated  over 
all  stations. 

5.  Station  detection  of  events  in  the  scattered 
population  is  simulated  by  first  computing  the  prob- 

A 

ability  of  detection  of  each  event  based  on  mb : 
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(3) 


k = 1,  2,  ... ,N(mb  ) 

2 1 
where  $(x)  = — i—  fXoc  e * ^dt 
/H 

m,  detection  threshold  (50%  level)  for  the 
b 

station  under  consideration. 

signal  to  noise  ratio  required  for  detection 
(same  for  all  stations). 

standard  deviation  of  the  noise  rms  for  mb 
for  the  station  under  consideration. 

Then  for  each  of  the  Pb  , a random  number  £k  uniformly 
distributed  between  0 and  1 is  computed,  and  the  two 
are  compared.  If  > Ph  the  event  is  considered 
"not  detected",  and  in  network  type  I the  (mb,  Mg) 
pair  is  eliminated  from  the  population.  If,  in  network 
type  I,  the  station  fails  to  detect  any  events  on  mb , 
a message  is  printed  and  that  station  is  omitted  from 
the  computation  of  network  detection  ctatistics. 

In  Network  type  I a station  will  not  detect  Mg 
unless  it  has  detected  mb>  This  might  correspond  to 
a network  made  up  from  the  bulletins  of  10  LASA's. 
Network  type  II  is  more  common  in  which  Mg  will  be 
looked  for  at  each  station  if  there  is  a network  detec- 

A * 

tion  on  mb*  In  network  type  II  we  retain  the  (mb,  Mg) 
pair  for  further  analysis  even  if  there  is  no  detection 


m 


bo 


R = 


n, 
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on  m.  . Also  at  this  point,  in  preparation  for  network 
l>  t h 

detection  statistics,  if  the  J event  is  detected  on 

the  element  NOMBD  (J)  is  increased  by  1 and  the 

array  AMBM  (J)  is  increased  by  m^ . 

After  this  "weeding”  process  is  completed,  this 
new  population  of  "observed"  events  is  searched  as  in 
step  3 to  determine  the  number  of  events  remaining  in 
each  magnitude  interval.  This  search  is  based  once  on 
the  M.  values  of  the  events,  and  again  on  the  in^  values. 
The  data  are  printed,  stored  individually,  and  also 
accumulated  over  all  stations. 

6.  For  each  station  the  events  observed  by  detection 
on  m,  at  that  station  (or  all  events  in  the  case  of 

b a 0 

network  type  II)  are  subjected  to  detection  on  M by 
the  same  method  used  in  step  5.  After  this  second 
"weeding"  process,  the  remaining  population  is  again 
searched  for  the  number  of  events  left  in  each  magni- 
tude interval.  These  data  arc  printed,  stored  indivi- 
dually and  accumulated  over  all  stations.  Again,  in 
preparation  for  netwrok  statistics  the  arrays  NOMSO(J) 

A 

and  AMSM(J)  are  increased  by  1 and  M respectively  if 

A 

detection  occurs  on  M . 

7.  Leas t- squares  regression  lines  are  fitted  to  the 
doubly  detected  events,  cither  to  the  entire  remaining 
population,  or  to  a truncated  portion  defined  by: 

Mil  < mb  < M2 2 MS  1 1 < Mg  < MS22 

where  the  limiting  values  are  specified  input  parameters. 
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Regression  of  both  Mg  on  m^  and  on  Mg  are  performed; 
moreover,  a maximum  likelihood  linear  fit  is  computed 
for  Mg-m^  as  in  step  2.  The  slopes  and  intercepts 
are  printed  and  accumulated  for  network  averages. 

8.  Step  4 is  repeated  for  the  observed  population. 

9.  The  ratio  of  the  number  of  events  detected  on 

both  m,  and  to  the  number  detected  on  m,  alone  is 

ns  b 

calculated  and  printed  for  each  magnitude  interval, 
first  using  M as  the  magnitude  of  each  event,  and 
again  using  m^.  These  ratios  represent  the  probability 
of  detection  on  given  prior  detection  on  m^(according 
to  type  I or  type  II  network)  as  a function  of  M and 

A ^ 

m^  respectively  at  the  given  station. 

10.  Finally,  a printer  plot  of  Mg-m^  for  the  events 

detected  on  both  m^  and  Mg  is  output  for  the  station. 

This  plot  also  shows  the  true  M -m,  line. 

s b 

Steps  1-10  are  repeated  for  each  station  in  the 
network.  After  all  stations  have  been  simulated,  the 
accumulated  data  is  presented  in  a network  summary 
for  the  first  set  of  events: 

Station  Summary 

I.  Over  all  initial  "scattered"  populations: 

1.  Average  slope  and  intercept  of  regressions 

A 

on  m^ . 

2.  Average  slope  and  intercept  of  regressions 

on  M . 

s 
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3.  Average  slope  and  intercept  of  maximum  like- 
lihood fits  of  M -m,  . 

s D 

4.  Accumulated  number  of  events  in  each  magni- 
tude interval  of  , both  incremental  and 
cumulative . 

5.  Accumulated  number  of  events  in  each  magni- 
tude interval  of  m^ , both  incremental  and 
cumulative. 

6.  Average  ni,  for  fixed  M . 

d s 

7 » Average  A for  fixed  m,  . 

s b 

II.  After  detection  on  m,  : 

b 

8.  Accumulated  number  of  events  in  each  magnitude 
interval  of  (same  as  1-4  in  the  case  of 
network  type  II). 

9.  Accumulated  number  of  events  in  each  magnitude 
interval  of  m,  . 

D 


After  detection  on  both  m,  and  M 

D S 

A 

(or  on  M only 

in 

the  case  of  network  type  II): 

10. 

Accumulated  number  of  events 

in  each  magnitude 

interval  of  M . 

s 

11. 

Accumulated  number  of  events 

A 

in  each  magnitude 

interval  of  m,  . 

D 

12. 

Average  m,  for  fixed  M . 

D s 

13. 

Average  M for  fixed  m,  . 

3 D 
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14. 


The  probability  of  detection  on  M given 

a 

prior  detection  on  m^  (or  without  considera- 
tion of  detection  on  m^  for  network  type  II) 
is  obtained  as  a function  of  Mg  by  computing 
the  ratio  of  the  data  of  step  10  to  the  data 
of  step  8 for  each  magnitude  interval. 

15.  The  probability  of  detection  on  M given 
prior  detection  on  m^  (same  as  abovej  is 
obtained  as  a function  of  m^  by  computing 
the  ratio  of  the  data  of  step  11  to  the  data 
of  step  9 for  each  magnitude  interval. 

16.  Average  over  stations  of  slope  and  intercept 

A A 

of  regression  fits  of  Mg  on  m^. 

17.  Average  over  stations  of  slope  and  intercept 

A A 

of  regression  fits  of  m^  on  Ms» 

18.  Average  over  stations  of  slope  and  intercept 

A A 

of  maximum  likelihood  fits  of  M and  n,  . 

s b 

Network  Summary 

A 

19.  The  probability  of  network  detection  on  m^ 
by  NRMBD  or  more  stations  as  a function  of 

m,  , is  calculated  by  seeing  if  N0MBD( J)>NRMBD 
(J  n event)  and  if  so  incrementing  by  1 the 
NMBNTD  cell  determined  by  SQGB=AMBM(J) /NOMBD 
(J)  = mb.  If  NOMBD (K)  < NkMBD  go  to  next 
event.  After  going  through  all  events, 

NMBNTD (K) /N  (m^)  is  printed  out. 
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20.  For  each  event,  after  (19)  if  NOMBD(J) >NRMBD , 
if  N0MSB_>NRMSD  then  we  increment  by  1 

the  NMSNTD  (K)  determined  by  SQGS=AMSM(K)  / NOMSD 
(K)  « Mg  and  the  NMBMSI)  (K)  cell  determined 
by  SQGB.  After  going  through  all  events 
NMSNTD  (K)/N  (m^)  is  printed  out.  Note  that 
because  of  magnitude  scatter  the  quantities 
printed  out  in  (19)  and  (20)  may  be  > 1.0. 

21.  To  obtain  the  probability  of  detection  of 
LR  by  at  least  NRMSD  stations  as  a function 
of  mb  we  compute  and  print  out  NMBMSD  (K)/ 
NMBNTD  (K).  These  numbers  are  always  < 1.0. 

22.  NMBNTD  (K)  , NMSNTD (K)  , and  NMBMSD  (K)  are 
printed  out,  together  with  their  cumulative 
versions. 

23.  Max-like  and  m^  and  M regression  lines  are 
computed  on  the  surviving  (mb,  ^g)  network 
pairs  for  NSTAT  (up  to  10)  sets  of  Mg  and 
mb  magnitude  intervals  (NSTAT)  x (EMI , EM 2 ; 

EMS 1 , EMS2) . 

24.  The  surviving  mb , $s  are  printed  out  and 
plotted  on  the  printer. 

25.  The  steps  19-24  are  repeated  for  NRDT  (up 
to  10)  pairs  of  (NRMBD , NRMSD)  values. 
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Final  Summary 

26.  The  entire  program  loops  NREP  times  (K=l, 
NREP);  saving  and  accumulating  the  max-like 
slopes  and  intercepts  as  well  as  NMBNTD  (K)  , 
NMSNTD  (K),  and  NMBMSD  (K)  in  KMBNTD  (K,  L) 
KMSNTD  (K,  L),  and  KMBMSD  (K,  L)(L=1,  NRDT) , 
When  all  repetitions  are  completed,  average 
regression  lines  are  printed,  along  with  the 
ratio  KMBMSD  (K,  L) /KMBNTD  (K,  L) . KMBNTD 
(K,  L),  KMSNTD  (K,  L) , and  KMBMSD  (K,  L)  are 
also  printed,  both  incremental  and  cumulative. 
Finally,  NRDT  printer  plots  showing  the 
accumulation  of  all  previous  network  plots 
for  each  L are  computed  and  printed. 

Correlated  Noise 

27.  If  NCORR  / 0 a preliminary  run  of  program 
SIMCRN  is  required  to  generate  noise  values 
with  the  required  correlation  matrix.  These 
are  output  to  a tape  which  must  be  mounted 
on  SYS1  when  run  on  the  360/44.  The  program 
is  iterated  as  many  times  (NREP)  as  it  is 
desired  to  iterate  PROGRAM  MSBNET.  Input 
for  program  SIMCRN  is  described  at  the  end 
of  this  writeup. 

4.  INPUT : Six  general  data  cards  followed  by  NSTAT 

cards  for  different  data  sets  for  regression  cal- 
culations and  by  one  data  card  for  each  station 
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(max.  30  stations)  are  required  for  each  model 
run  of  program  MSBNET . The  variables  and  formats 
of  the  cards  are  as  follows: 

CARD  1 (6110,  415) 

IX  Initial  value  for  random  number  genera- 

tion. (If  zero,  IX  is  set  to  777773777.) 

I PNT  If  zero,  omit  printout  of  individual  sta- 

tion data.  Network  summaries  and  the 
final  summary  will  still  be  printed. 

IHPKL  If  zero,  omit  printout  of  individual  Mg, 
event  pairs  for  each  station. 

JIPKL  If  zero,  omit  printout  of  Mg , m^  event 
pairs  for  network. 

JZ3Z  If  zero,  network  type  I.  If  cne,  network 
type  II. 

NRDT  Number  of  (NRMBD.NRMSD)  pairs  to  be  used 
(see  cards  3 and  4) . 

MBMIN  Limits  on  the  and  Mg  axes  of  the 

MBMAX  printer  plots  of  station  and  network  data. 

MSMIN  Default  to  2. 0-7,0  on  both  axes.  If 
MSMAX  MBMIN  * -1,  the  maximum  and  minimum  data 
values  will  be  used. 

CARD  2 (6F10.2,  2110) 

Ml  Smallest  "true"  m,  value  used  in  the 

b 

generation  of  each  population  of  "observed" 
events . 
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A Parameters  of  the  seismicity  density  func- 

ALP  tion  (see  cq.  1)  must  be  adjusted  to  give 

no  more  than  5000  events:  A/ALP<5000. 

BBT A Slope  and  intercept  of  the  "true"  Mg-mb 

Cl  NT  relation:  Mg  = BETA  * m^  + CIN’T. 

R Signal  to  noise  ratio  required  for  detection. 

N Number  of  stations  (must  not  exceed  30). 

NREP  Number  of  repetitions. 

CARD  3 (1015) 

NRMBD(I)  Minimum  number  of  station  detections  of 
mb  required  to  declare  a network  detec- 
tion. (1=1,  NRDT ; see  card  1) 


CARD  4 (1015) 

NRMSD(I)  Minimum  number  of  station  detections  of 
M§  required  to  declare  a network  dctec- 
tion0  (I  = 1,  NRDT;  see  card  1) 

CARD  5 (4F10.4,  110) 

SGNTB  Standard  deviation  of  network  mb  for  maxi- 
mum likelihood  regressions. 

SGNTS  Standard  deviation  of  network  M for  maxi- 
mum likelihood  regressions. 

SIGSS  Standard  deviation  for  mb  source  bias. 

SIGBS  Standard  deviation  for  M source  bias. 

NCORR  Switch  0)  for  a correlated  noise  and 
signal  run. 


- AI 1 1 - 1 3- 


CARD  6 (4F10.2,  110) 


Mil  Lowest  magnitude  mb  to  consider  for  regres- 
sion fit  of  detected  events. 

Mi,2  Highest  magnitude  m^  to  consider  for  regres- 
sion fit  of  detected  events. 

MS 11  Lowest  magnitude  M^  to  consider  for  regres- 
sion fit  of  detected  events. 

MS 2 2 Highest  magnitude  Mg  to  consider  for  regres- 

sion fit  of  detected  events. 

NSTAT  Number  of  sets  of  EMI,  EM2 , EMC1,  EMS2  to 
follow.  If  NSTAT  = 0,  station  data  cards 
come  next. 


CARD  7 - CARD  (6  ♦ NSTAT)  (4F10.2) 

EM1  C1)  Truncation  limits  on  mb 

EM2  (I)  j _ ^ NSTAT  anc*  Ms  ^or  ac*ditional  re- 
EMS1(I)  gression  fit  (same  defini- 

LMS 2(1)  tions  respectively  as 

Card  6) . 


Ith  STATION  DATA  CARD  (I  = 1,  N)  (5F10.4) 

AMBO(I)  mb  detection  threshold  for  the  Ith  station, 
50%. 

AMSO(I)  M detection  threshold  for  the  Ith  station 
50%. 


ASGMB(I)  Standard  deviation  of  observed  events 
about  the  Mg-mb  relation  along  the  mb 
axis  for  the  Jth  station,  with  variance 
due  to  source  bias  excluded. 
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ASGMS(I)  Standard  dcviat  ion  of  observed  events 

about  the  M.-m,  relation  alone  the  M 
s d s 

axis  for  the  Ith  station  with  variance 
due  to  source  bias  excluded. 

ASGB(I)  Standard  deviation  of  the  noise  rms  for 

f h 

mb  at  the  I station. 

ASGS(I)  Standard  deviation  of  the  noise  rms  for 
Mg  at  the  station. 

Iirror  Messages: 

Self-explanatory. 

Memory  requirements: 

200  K on  360/44.  Makes  extensive  use  of  the 
auxiliary  disc  memory. 

References:  Seismic  Threshold  Determination,  TR-74-3, 

D.  H.  von  Seggern,  R.  R.  Blandford,  1974. 

Programmers:  11.  Dusted,  T.  Mchlfrcsh,  R.  Blandford. 
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C PROGRAM  SIM  CRN  CORRELATED  DATA  FOR  NETWORK 
C PROGRAM  TO  GENERATE  5000  CORRELATED  OBSERVATIONS  FOR  10  STATIC 
DIMENSION  C0RR(  10,5000) ,CCVAR(  10 , 10) ,AMEAN( 10) ,X( 10) 
DIMENSION  Y( 10) 


I- 5  I FC0 R=  1 COVARIANCE  MATRIX  IS  INPUT 

* 2 CORRELATION  MATRIX  IS  INPUT 

6-10  N NUMBER  OF  STATIONS  USED  IN  SIMULATION  (N) 

II- 20  N RAN  STARTING  RANDOM  NUMBER  (DEFAULT  TO  77773777) 
21-30  NREP  NO.  OF  TIMES  TO  REPEAT  SIMULATION 

1-10  AMEAN(l)  THE  MEAN  OF  THE  FIRST  STATION 

71-80  AMEAN( 8)  THE  MEAN  OF  THE  8TH  STATION 

REPEAT  UNTIL  ALL  STATIONS  MEANS  ARE  READ 

1-5  STANDARD  DEVIATION  OF  STATION  1 


76-80  STANDARD  DEVIATION  OF  STATION  16 

REPEAT  UNTIL  ALL  STATION  STANDARD  DEVIATIONS  ARE  READ 


ONLY  INPUT  THE  LOWER  TRIANGULAR  PORTION  OF  MATRIX 
CARD  1 AS ( 1 , 1 ) 

CARD  2 AS ( 2 , 1 ) , AS( 2 ,2 ) 

CARD  3 AS (3.1)  , AS ( 3 , 2 ) , AS ( 3 , 3 ) 


C 

C 


CARD  NST  AS ( NST  , 1 ) ,AS(NST,2) , AS ( NS  T , 3 ) 


.AS(NST.NST) 
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March  1,  1973 


A SUBROUTINE  FOR  MAXIMUM  LIKELIHOOD  LINEAR  FITTING 

ON  THE  IBM  360-44 


1.  Subroutine  Name:  MAXLIK 

2.  Purpose : 

This  subroutine  performs  maximum  likelihood  fitting 
of  a straight  line  to  measured  points  having  correlated 
and  normally  distributed  errors  in  both  variables.  The 
method  involves  point  adjustments  along  likelihood 
ellipse  diameters  conjugate  to  the  fitted  line. 


3.  Calling  Sequence: 

CALL  MAXLIK  (X,  Y,  N,  G,  R,  A,  B,  SA , SB,  SX,  SY) 

4 . Input : 


X = 
Y = 
N = 
G = 
R = 


1st  variable  array 
2nd  variable  array 
Number  of  data  points  in  X, 
Ratio  of  standard  deviations 
Correlation  coefficient  of  e 


Y arrays 

a /a 
y'  x 

rrors  in  X and  Y 


5 . Output : 

A => Slope  of  the  fitted  line:  Y = A*  X + B. 

B = Intercept  of  the  fitted  line:  Y = A*  X + B. 
SA  = Standard  deviation  of  the  slope  for  the 

uncorrelated  (R  = 0.)  case.  If  R + 0,  the 
returned  value  of  SA  is  -1.0. 
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SB  = Standard  deviation  of  the  intercept  for  the 
uncorrelated  (R  = 0.)  case.  I f R i 0 , the 
returned  value  of  SB  is  -1.0. 

SX  = Standard  deviation  of  the  X observations  with 
respect  to  their  true  values. 

SY  = Standard  deviation  of  the  Y observations  with 
respect  to  their  true  values. 

Error  Messages:  None. 

Memory  Requirements:  Less  than  2K  bytes 

References : 

Ericsson,  U.,  Maximum  likelihood  linear  fitting 

when  both  variables  have  normal  and  correlated 
errors;  The  Research  Institute  of  National 
Defense,  Stockholm  80,  Sweden,  FOA  4 Rapport, 
C4474-A1 , Sept.  *71. 
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MAXLIK  DATE  03/02/73 

SUBROUTINE  MAX L I K( X , Y ,N ,G , RHC , A,B ,S A ,SB ,S X ,S Y ) 

C 

DIMENSION  X( 1)  ,Y( 1) 

C 

TN=FLOAT ( N) 

C 

XB=0.0 
YB=0.0 
DO  30  J = 1 , N 
XB  = XB  + X ( J ) 

YB  = YB  + Y ( J ) 

30  CONTINUE 
XB^XB/TN 
YB=YB/TN 
C 

SUU=0.0 

SUV*0.0 

SVV=0.0 

DO  40  J-l.N 

U J = X ( J ) - X B 

VJ - Y ( J ) - YB 

SUU=SUU+(UJ*UJ) 

SUV=SUV+(UJ*VJ) 

SVV=SVV+( VJ*VJ) 

40  CONTINUE 
C 

BX=SUV/SUU 
AX  = XB- ( B X * XB  ) 

BY=SUV/S VV 
AY  = XB - ( B Y*YB  ) 

FAC-SQRT( 1. /(TN-2. ) )*( ( l./(BX*BY) )- 1 . ) ) 

jBX=BX*FAC 

SBY  = B Y *F  AC 

G2=G*G 

S0= ( RHO*B  Y ) - G2 
Sl  = ( G2/B X) - 1 . /BY  ) 

S2=1.-RH0*G*(  l./BX) 

C 

A=(-S1+SQRT(S1*S1-4.*S0*S2))/(2.*S2) 

B=YB-A*XB 

C 

3S=0.0 

DO  50  J = 1 , N 

F = Y ( J ) -B - A*X ( J ) 

SS=SS+( F*F) 

50  CONTINUE 

SX«SQRT(SS*(  I./(G2  + A*A-2.*A*RHO*G))*(l./(TN-2.))) 
SY=G*SX 
C 

IF(RHO.EQ.O.O)  GO  TO  60 

S A=-  1 . 0 

SB*-1.0 

RETURN 

60  DA»(G2+(BX/BY))/(TN-1.) 

DB=G2+A*A 

DC= ( 1 . / ( BX*B Y ) ) - 1 . 0 
SA=A*SQRT(DA*DB*DC) 

SB=SQRT((SVV/(TN*(TN-2.)))+(XB*XB*SA*SA)) 


C 


RETURN 

END 
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